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Abstract 

Modern  space  situational  awareness  is  focused  on  the  detection,  tracking,  identification,  and  char¬ 
acterization  of  passive  and  active  resident  space  objects.  In  the  past,  this  process  relied  primarily  on 
ground-based  sensors.  However,  difficulties  arise  when  smaller,  more  distant,  or  otherwise  dim  objects 
are  considered.  Space  based  sensors  may  help  alleviate  some  of  these  difficulties.  The  work  funded  under 
this  award  has  focused  primarily  on  optimal  design  for  a  constellation  of  space  based  sensors  to  support 
the  goals  of  space  situational  awareness.  To  that  end,  the  work  leverages  concepts  from  computer  graph¬ 
ics  to  analyze  sensor  coverage  strictly  from  a  numerical  perspective,  enabling  a  unique  design  approach 
that  is  adaptable  to  any  sensor  network  for  situational  awareness  applications. 


1  Executive  Summary 

As  near-Eartli  space  becomes  increasingly  crowded  with  spacecraft  and  debris,  the  need  for  improved  space 
situational  awareness  (SSA)  has  become  paramount.  Modern  SSA  is  concerned  with  the  detection,  tracking, 
identification,  and  characterization  (DTI&C)  of  passive  (e.g.  debris)  and  active  (i.e.  maneuverable)  resident 
space  objects  (RSOs),  at  all  altitudes,  with  known  accuracy  and  precision.  Modern  advances  in  technology 
miniaturization,  and  the  growing  capabilities  of  nano-  and  pico-satellite  platforms,  have  raised  many  new  and 
exciting  technical  challenges  in  this  area.  Contemporary  ground-based  systems,  such  as  the  space  surveillance 
network  (SSN),  are  oftentimes  unable  to  resolve  small,  distant,  or  otherwise  faint  objects.  In  these  instances, 
space-based  sensors  may  be  beneficial  to  augment  existing  SSA  capabilities  for  enhanced  DTI&C. 

In  this  study  SSA  goals  are  primarily  accomplished  through  a  constellation  of  space  based  sensors,  though 
the  approach  also  allows  for  the  inclusion  of  ground  based  sensors.3  The  constellation  is  designed  to  provide 
optimal  above-the-horizon  (ATH)  coverage.  This  is  fundamentally  different  from  more  classical  optimal 
constellation  coverage  problems,  which  focus  on  optimizing  coverage  of  ground-based  targets.  The  latter  is 
typically  referred  to  as  the  below-the-horizon  (BTH)  coverage  problem.  The  ‘horizon,’  in  this  case,  is  defined 
as  a  line  drawn  from  the  satellite  and  tangent  to  the  surface  of  a  fictitious  sphere,  one  concentric  with  the 
Earth  but  of  equal  or  greater  radius.  This  reference  tangent  height  sphere  is  simply  a  convenient  construct 
that  may  be  used,  for  instance,  to  represent  the  bounds  of  the  atmosphere.  From  a  three-dimensional 
perspective,  this  line  traces  a  cone  that  is  tangent  to  this  reference  sphere.  RSO’s  that  exist  above  this  cone 
are  said  to  be  viewed  against  a  space  background.  Thus,  a  sensor  concerned  with  such  targets  is  said  to 
provide  ATH  coverage. 

Aside  from  being  primarily  interested  in  RSO’s  that  exist  above  the  horizon,  the  present  study  is  further 
interested  in  objects  that  evolve  within  some  prescribed  three-dimensional  region  of  space  that  exists  above 
the  horizon.  The  goal  of  this  effort  is  to  identify  a  generalized  numerical  process  by  which  to  accomplish  the 
optimal  design  of  a  constellation  of  orbiting  sensors  that  meet  a  specific  collective  objective  for  ATH  cov¬ 
erage.  This  includes  consideration  of  non-coplanar  elliptical  orbits,  arbitrary  sensor  region’s  of  regard,  and 
dynamical  effects.  The  modeling  approach  selected  considers  each  sensor’s  region  of  regard  and  any  imposed 
visibility  constraints  as  time  varying  and  dynamic  surfaces.  Thus,  the  coverage  provided  by  the  constellation 
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can  be  thought  of  as  a  time- varying  volume,  which  is  subsequently  defined  as  the  cost-index  for  an  optimiza¬ 
tion  process.  The  calculation  of  this  time- varying  volume  of  coverage  is  accomplished  by  performing  a  series 
of  boolean  operations  between  the  reference  surfaces  at  each  time  step  during  the  simulation. 

As  an  initial  step,  the  first  year  of  this  award  focused  on  investigating  the  preliminary  aspects  of  the 
simpler  planar  constellation  design  problem.  A  two-pronged  approach  was  selected,  one  focusing  on  a 
strictly  numerical  approach,  and  another  focusing  on  devising  a  closed  form  expression  for  constellation 
coverage  that  could  be  used  in  an  optimal  constellation  design  process.  Investigating  possible  numerical 
techniques  initially  revealed  polygon  clipping  algorithms  as  a  suitable  starting  point.  Relevant  numerical 
methods  were  surveyed,  studied,  evaluated,  implemented,  and  validated  against  publicly  available  analytical 
examples  for  the  single  satellite  ATH  coverage  problem.  The  analytical  approach  focused  on  investigating 
the  mathematical  preliminaries  involved  in  a  purely  geometric  approach. 

In  the  second  year  of  this  effort,  the  ATH  coverage  provided  by  planar  constellations  was  investigated 
in  depth  from  both  the  numerical  and  analytical  perspectives.  The  analytical  approach  focused  on  the 
mathematical  aspects  involved  in  addressing  coverage  multiplicities.1,2  The  numerical  approach  focused  on 
the  development  of  a  numerical  process,  based  on  polygon  clipping  techniques,  to  numerically  analyze  planar 
constellations  for  ATH  coverage  of  any  desired  multiplicity  (i.e.  number  of  satellites  that  can  observe  a  given 
region  simultaneously).3  The  numerical  approach  uses  sequences  of  Boolean  operations  to  determine  the 
multiple  overlap  of  polygons  representing  the  in-plane  sensor  coverage  regions  of  each  satellite  and  the  target 
region.  The  advantage  of  restricting  the  investigation,  at  least  initially,  to  planar  constellations  of  satellites 
sensing  in  regions  symmetric  across  the  orbit  plane  is  that  the  in-plane  ATH  coverage  region  is  proportional 
to  volumetric  ATH  coverage.  Thus,  by  addressing  the  2D  problem  one  indirectly  addresses  the  3D  problem 
as  well.  The  numerical  ATH  coverage  algorithm  was  applied  to  several  example  design  constellation  problems 
using  a  Non-Linear  Programming  (NLP)  solver  to  illustrate  how  it  may  be  used  in  a  constellation  design 
process. 

The  final  year  of  this  award  initially  focused  on  closing  out  the  planar  analysis  to  begin  the  transition  to 
the  full  three-dimensional  problem.  During  that  time,  the  analytical  constellation  analysis  methods2  devised 
were  successfully  employed  in  validating  the  numerical  process.  To  further  demonstrate  the  integration 
of  the  numerical  coverage  calculation  with  an  on-line  optimization  process,  a  Mixed  Integer  Non-Linear 
Programming  (MINLP)  algorithm  was  applied  to  several  example  constellation  design  problems,  expanding 
the  types  of  planar  constellation  optimization  problems  that  may  be  addressed.4 

As  implemented,  extensions  to  the  three-dimensional  case  entail  discretizing  volumes  representing  cover¬ 
age  and  target  regions  into  a  prescribed  number  of  planar  cross-sections  of  uniform  separation.  The  coverage 
area  within  each  of  the  planar  cross-sections  is  evaluated  using  the  previously  developed  planar  analysis 
algorithm,  and  is  subsequently  multiplied  by  the  cross-section  plane  separation  distance.  This  approximates 
the  volume  in  each  ‘layer,’  that  is  summed  with  the  volume  of  every  other  layer  of  a  desired  coverage  multi¬ 
plicity  to  approximate  the  coverage  volume.  This  necessitates  a  method  to  generate  the  cross-sections  of  the 
coverage  volumes,  which  is  straightforward  for  the  omni-directional  sensor  profile  case.  This  new  approach 
is  applied  to  several  example  problems  demonstrating  volumetric  ATH  coverage  evaluation  of  constellations 
of  non-coplanar  omni-directionally  sensing  satellites.  Since  the  results  of  all  related  initial  investigations, 
during  the  first  two  years  of  this  award,  have  already  been  published  in  conference  articles,  graduate  theses, 
or  are  pending  submission  to  specific  archival  journals,  the  present  report  offers  only  a  relevant  summary  of 
those  activities  but  otherwise  focuses,  in  later  sections,  to  the  details  of  the  three-dimensional  extension. 

2  Accomplishments/New  Findings 

All  optimization  processes  require  that  a  cost  index  or  performance  metric  be  pre-defined.  Typically,  such 
metrics  are  defined  as  exact  or  closed  form  expressions,  an  element  that  is  certainly  required  in  any  numerical 
solution  process.  That  is,  perhaps,  the  most  difficult  aspect  of  the  problem  investigated  during  the  course 
of  this  study.  The  results  of  this  work  led  to  the  design  of  a  numerical  process  that  allows  for  the  coverage 
volume  provided  by  a  constellation  to  be  systematically  computed  numerically.  The  process  treats  the 
region  of  regard  of  each  satellite,  the  region  of  interest,  and  the  visibility  constraints  as  a  series  of  surfaces. 
The  coverage  volume  is  identified  through  a  series  of  boolean  operations  between  the  reference  surfaces  in 
the  constellation.  The  methodology  is  integrated  with  various  optimization  methods  to  demonstrate  the 
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integrated  process  by  which  the  optimal  orbital  design  is  accomplished.  The  innovation,  in  this  case,  is  not 
in  the  optimization  process  itself.  Rather,  it  is  in  generalized  numerical  process  designed  to  address  the 
calculation  of  the  cost  index. 

The  results  of  this  work  can  be  applied  in  many  contexts,  from  constellation  design  to  close-proximity 
operations.  Thus,  far,  the  results  of  this  work  have  been  published  in  two  conference  papers,  two  MS  theses. 
Two  journal  articles  are  currently  in  preparation  for  submission  to  the  AIAA  Journal  of  Spacecraft  and 
Rockets.  The  grant  also  partially  funded  the  initial  year  of  a  follow-on  PhD  level  research.  The  publications 
previously  listed  are  summarized  below: 

•  A.T.  Takano  and  B.G.  Marchand,  ’’Optimal  Constellation  Design  for  Space  Based  Situational  Aware¬ 
ness  Applications,”  AAS/AIAA  Astrodynamics  Specialists  Conference,  Girdwood,  AK,  August  2011. 
Paper  No.  AAS11-543 

•  A.D.  Biria  and  B.G.  Marchand,  ’’Constellation  Design  for  Space-Based  Situational  Awareness  Applica¬ 
tions:  An  Analytical  Approach,”  2011  AAS/AIAA  Astrodynamics  Specialists  Conference,  Girdwood, 
AK,  August,  2011.  Paper  No.  AAS11-538. 

•  Numerical  Analysis  and  Design  of  Satellite  Constellations  for  Above  the  Horizon  Coverage,  Andrew 
Takano,  MS  Thesis,  Advisor:  Dr.  Belinda  Marchand,  December  2010. 

•  Analytical  Approach  to  the  Design  of  Optimal  Satellite  Constellations  for  Space-Based  Space  Situa¬ 
tional  Awareness,  MS  Thesis,  Advisor:  Dr. Belinda  Marchand,  December  2011. 

The  numerical  methodology  devised  has  been  validated  with  an  analytical  approach,1,2  also  devised 
during  this  investigation  under  a  simplified  set  of  assumptions.  The  analytical  approach,  however,  is  only 
applicable  to  the  case  of  planar  constellations,  where  satellites  are  collocated  along  the  same  circular  orbit. 
The  level  of  complexity  associated  with  the  general  three-dimensional  case,  however,  is  best  addressed  by  the 
numerical  process  devised  during  this  study.  The  numerical  approach  has  been  successfully  demonstrated 
in  the  solution  of  a  MINLP  optimization  problem  in  the  planar  case,  where  the  goal  was  to  determine  the 
optimal  number  of  satellites  in  the  constellation.  Also,  the  numerical  process  has  been  demonstrated  in  a 
volumetric  coverage  example  for  non-coplanar  elliptical  orbits  of  non-collocated  satellites.  At  the  time  of  the 
conclusion  of  this  grant,  only  omni-directional  sensors  could  be  considered  in  the  full  three-dimensional  case. 
However,  non-omnidirectional  sensors  were  successfully  addressed  in  the  two-dimensional  case  during  the  first 
two  years  of  this  award.  Extending  the  generalized  numerical  process  to  accommodate  non-omnidirectional 
sensors  in  the  three-dimensional  case,  as  well  as  relative  attitude,  and  other  complexities,  is  a  subject  that 
requires  further  study.  Another  aspect  that  requires  further  study  is  the  numerical  efficiency  of  the  overall 
process.  A  parallelized  structure  that  leverages  GPU’s  is  envisioned  as  a  possible  follow-up  path  that  would 
allow  for  a  more  computationally  efficient  process  that  could  effectively  address  these  added  complexities. 
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3  Project  Background 

The  problem  of  interest  in  this  investigation  pertains  to  optimal  sensor  coverage  of  targets  against  a 
space  background  as  provided  by  a  generalized  network  of  orbiting  sensor  platforms.  As  a  whole,  the 
constellation  offers  a  dynamic  heterogeneous  network  of  sensors,  each  contributing  a  unique  set  of  capabilities. 
A  candidate  optimal  design  process  may  seek  to  optimize  over  a  set  of  sensor  specific  parameters,  a  set  of 
orbital  parameters,  or  both.  In  this  work,  the  cost  metric  employed  is  the  overall  coverage  provided  by 
the  constellation  within  the  region  of  interest.  In  the  modeling  approach  selected,  the  sensor  coverage,  the 
region  of  interest,  and  all  related  visibility  constraints  are  modeled  as  dynamically  evolving  three-dimensional 
surfaces. 

Earlier  studies,  by  Marcliand  and  Kobel,5  considered  a  simplified  version  of  this  problem.  Specifically, 
for  a  single  satellite,  their  goal  is  to  maximize  sensor  visibility  of  targets  against  a  space  background  when 
the  region  of  interest  is  bounded  between  an  upper  and  lower  target  altitude.  The  results  of  this  study  are 
critical  to  the  present  investigation  because  they  establish  a  meaningful  and  easily  verifiable  baseline  for 
comparison. 

The  goal  of  this  effort  was  to  devise  a  generalized  numerical  approach  that  would  enable  the  optimal 
design  and  analysis  of  constellations  aimed  at  providing  some  measure  of  above-the-horizon  coverage.  This 
process  should  allow  for  consideration  of  heterogeneous  sensors  arbitrarily  distributed  among  a  series  of 
non-coplanar  elliptical  orbits.  Naturally,  this  is  a  very  broad  specification  that  adds  an  immense  degree 
of  complexity  to  the  problem.  Attempting  to  identify  a  closed  form  solution  under  these  circumstances  is 
simply  not  feasible.  This  motivates  the  numerical  approach  proposed. 

In  implementing  the  numerical  approach,  it  is  assumed  that  the  region  of  regard  of  any  given  sensor 
in  the  network  can  be  represented  as  an  arbitrary  three-dimensional  sensor  volumetric  mesh  (SVM).  One 
advantage  of  this  approach  is  that  the  user  need  not  have  any  a  priori  knowledge  of  the  actual  sensor 
hardware.  Consider,  then,  that  each  satellite  in  a  constellation  contributes  at  least  one  SVM  to  the  network. 
Furthermore,  the  attitude  of  each  satellite  in  the  constellation  may  not  necessarily  be  fixed  in  relation  to 
the  reference  frame  of  interest.  Thus,  the  SVM  is  possibly  both  translating  and  rotating  in  relation  to  the 
region  of  interest.  Subsequently,  as  a  given  satellite  orbits  along  its  path,  its  SVM  may  or  may  not  overlap 
with  the  ATH  region  of  interest.  Determining  the  collective  coverage  provided  by  the  constellation  at  any 
given  instant  of  time,  then,  is  reduced  to  the  following  questions: 

•  Which  of  the  satellites  in  the  constellation  are  passing  through  or  are  near  the  region  of  interest  at 
that  time? 

•  Of  those  satellites,  how  many  SVM’s  are  properly  positioned  and  oriented  so  as  to  intersect  the  region 
of  interest? 

•  What  is  the  volume  of  the  intersection  of  each  SVM  with  the  region  of  interest? 

•  What  is  the  combined  volume  provided  by  the  union  of  these  SVM  intersections  with  the  region  of 
interest? 

•  What  is  the  overlap  volume  (if  any)  between  neighboring  SVM’s  within  the  region  of  interest? 

•  Is  the  entire  region  of  interest  covered  by  these  satellites?  or  are  there  gaps  in  coverage? 

•  If  gaps  are  present,  how  can  the  constellation  parameters,  or  the  attitude  of  the  vehicles,  be  adjusted 
such  that  coverage  is  maximized? 

These  questions  are  all  addressed  at  each  instant  of  time  during  a  simulation  or  optimization  process  since 
the  SVM’s  of  each  satellite  are  dynamically  evolving  in  position  and  attitude.  Optimizing  the  constellation, 
then,  may  involve  the  identification  of  orbital  parameters  that  maximize  the  coverage  of  that  area  over  time. 
Alternatively,  it  may  require  that  one  determine  the  minimum  number  of  satellites  required,  and  their  orbital 
characteristics,  to  achieve  the  coverage  requirements  of  that  region. 

To  demonstrate  the  above  process,  a  generalized  algorithm  for  numerically  estimating  the  coverage  pro¬ 
vided  by  a  constellation  was  formulated  and  validated.  Since  analytical  results  are  only  available  in  the 
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simplified  two-dimensional  case,  the  validation  was  restricted  to  those  examples  before  proceeding  to  in¬ 
corporate  three-dimensional  extensions  of  the  work.  The  validation  was  based  on  the  results  published  by 
Marchand  and  Kobel5  and  Marchand  and  Biria.2  These  earlier  studies  focused  on  identifying  analytical 
expressions,  based  on  geometrical  arguments,  for  the  coverage  area  as  a  function  of  altitude.  The  region 
of  interest  was  bounded  by  an  upper  and  lower  target  altitudes.  Strictly  to  introduce  the  general  notation 
employed  in  this  report,  Figure  1  illustrates  a  cross-sectional  view  of  the  single  satellite  dual-altitude  band 
ATH  coverage  problem.5  In  this  figure,  the  upper  and  lower  target  altitude  shells  are  labeled  UTAS  and 
LTAS,  respectively.  The  satellite  is  identified  as  point  S.  The  horizon  is  defined  relative  to  the  tangent  height 
shell,  labeled  THS.  Both  of  these  earlier  investigations  specifically  focused  on  omnidirectional  sensors.  This 
is  represented  in  Figure  1  as  a  circular  cross  section  centered  at  the  satellite.  The  intersections  of  the  tangent 
lines  and  the  sensor  range  are  all  labeled  on  the  figure  as  well.  The  regions  shaded  in  yellow  represent  the 
ATH  region  covered  by  the  satellite  at  its  current  altitude,  and  that  area  is  a  piece- wise  differentiable  function 
of  the  satellite  altitude.5  Of  course,  in  the  generalized  case,  where  the  sensors  are  not  omni-directional  and 
the  satellites  in  the  constellation  are  unevenly  distributed  among  multiple  non-coplanar  orbits,  an  analytical 
solution  is  not  available.  This  motivates  the  numerical  approach  that  is  the  focus  of  this  effort. 
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4  Planar  Constellation  Analysis 

During  the  initial  phase  of  this  effort,  a  generalized  numerical  approach  is  devised  for  the  calculation  of  ATH 
coverage  with  any  desired  level  of  coverage  multiplicity.  The  approach  described  involves  discretization 
(into  polygons)  of  in-plane  satellite  coverage  regions  and  regions  of  interest  (assumed  to  be  a  dual-altitude 
band  shell  bounded  above  and  below  by  prescribed  altitudes).  The  interactions  between  these  polygons  are 
analyzed  by  performing  different  sequences  of  polygon  clipping  operations  (the  computational  method  used 
to  perform  Boolean  operations  between  coplanar  polygons),  yielding  a  result  polygon  that  represents  the 
region  exhibiting  the  desired  coverage  multiplicity  bounded  by  the  region  of  interest.  When  the  enclosed 
area  of  this  result  polygon  is  computed,  it  provides  an  effective  measure  of  in-plane  ATH  coverage  for  a  given 
constellation  configuration. 

The  necessary  sequences  of  polygon  clipping  operations  for  this  analysis  are  first  developed  using  set 
notation  with  a  focus  on  efficiency  by  avoiding  unnecessary  polygon  clipping  operations.  These  efficient 
ATH  coverage  algorithms  are  then  implemented  using  several  different  polygon  clipping  implementations 
which  are  compared  by  their  performance.  The  error  introduced  by  approximating  curvilinear  regions  with 
finite-resolution  polygons  is  extensively  analyzed,  and  relationships  guiding  appropriate  polygon  resolution 
selection  for  a  desired  accuracy  are  developed. 

Several  example  problems  are  then  solved  that  illustrate  how  the  numerical  ATH  coverage  algorithms 
can  be  incorporated  into  a  satellite  constellation  design  problem.  Using  a  parameter  optimization  code,  the 
coverage  algorithms  are  first  used  as  objective  functions  to  maximize  ATH  coverage,  then  as  constraints 
requiring  a  target  ATH  coverage  amount  be  achieved.  Simple  though  they  may  be,  these  examples  demon¬ 
strate  the  utility  of  the  numerical  ATH  coverage  algorithms  in  the  design  of  constellation  configurations  for 
ATH  coverage. 

It  is  worth  emphasizing  the  generality  this  methodology  allows  in  the  analysis  of  ATH  coverage.  First, 
the  sensor  cross-sections  can  be  of  any  shape,  and  can  even  be  unique  to  each  satellite  in  a  constellation. 
Secondly,  the  analysis  can  be  used  in  both  time-invariant  and  time-varying  analyses  -  the  only  two  time- 
varying  parameters  necessary  are  the  locations  and  in-plane  attitudes  of  each  satellite  at  a  given  time.  Lastly, 
although  the  study  only  considers  the  dual-altitude  band  shell  as  the  region  of  interest,  it  is  possible  to  analyze 
coverage  in  regions  of  interest  that  are  of  any  desired  planar  form.  For  instance,  the  region  of  interest  could 
be  fixed  above  a  certain  location  on  the  central  body,  i.e.  to  analyze  for  ATH  coverage  specifically  above  a 
prescribed  geographic  region. 

A  basic  example  of  planar  sensor  region  and  region  of  interest  discretization  is  presented  in  Section 
4.1.  Following  that,  set  notation  expressions  defining  the  sequences  of  Boolean  operations  are  developed 
in  Sections  4.2,  4.3,  and  4.4.  These  methods  determine  regions  of  single,  double,  and  arbitrarily  defined 
multiplicity  of  overlap  among  the  sensor  regions  within  the  region  of  interest. 

4.1  Discretization  of  Sensor  Regions  and  Region  of  Interest 

Omni-directional  satellite  sensor  profiles  are  the  simplest  possible  scenario,  and  are  a  straightforward  assump¬ 
tion  to  implement  in  the  planar  case.  Each  satellite  is  assumed  to  have  sensor  visibility  within  a  prescribed 
radius.  However,  the  sensors  are  also  assumed  to  only  be  capable  of  observing  targets  that  are  against  a 
space  background  from  the  perspective  of  the  satellite.  Thus,  any  portion  of  the  sensor  region  in  the  direction 
of  the  Earth  must  be  omitted  from  consideration  (it  is  a  blind  spot).  The  resulting  effective  ATH  coverage 
region  ( RSe )  provided  by  a  single  satellite  is  shown  in  Figure  2a.  The  omni-directional  sensor  assumption 
is  not  necessary  to  analyze  problems  using  the  methods  developed  in  this  study.3  It  does,  however,  allow  for 
a  significantly  more  straightforward  presentation,  and  is  consequently  maintained  in  the  current  discussion. 

The  target  region,  or  region  of  interest  (AS),  in  the  current  discussion  is  assumed  to  be  a  dual-altitude 
band  region,  as  shown  in  Figure  2b.  Lower  and  upper  bounds  in  altitude  are  prescribed  that  define  the 
range  of  altitudes  where  ATH  coverage  is  of  interest.  As  in  the  case  of  the  assumed  omni-directional  satellite 
sensor  profiles,  this  assumption  is  not  due  to  a  limitation  of  the  methods  developed  in  this  study,  but  is 
selected  to  simplify  the  present  discussion.  The  region  of  interest  can  be  of  any  planar  form,  although  the 
time  invariance  of  the  problem  is  lost  in  the  general  case. 

Without  loss  of  generality,  the  assumed  forms  of  RSe  and  AS  shown  in  Figure  2  are  used  to  illustrate 
the  coverage  expressions  developed  in  Sections  4.2,  4.3,  and  4.4. 


6 


(a)  Effective  ATH  Satellite  Sensor  Region  ( RSe ) 


(b)  Region  of  Interest  (AS) 


Figure  2:  Sensor  Region  and  Region  of  Interest 


4.2  Single  Coverage 

Once  the  individual  satellite  ATH  coverage  regions  and  region  of  interest  are  discretized  into  polygons,  the 
total  region  of  single  coverage  within  the  region  of  interest  can  be  determined  using  n  Boolean  operations 
(performed  using  polygon  clipping).  This  process  is  illustrated  in  Figure  3.  Figures  3a-3i  show  how  the  total 
effective  ATH  coverage  region  is  formed  using  n  —  1  union  operations.  Subsequently,  the  total  effective  ATH 
coverage  region  is  intersected  with  the  region  of  interest,  as  shown  in  Figures  3j-31,  yielding  the  region  of  ATH 
coverage  in  the  region  of  interest,  Cix.  The  area  of  this  polygon  region  is  readily  determined  numerically,6 
and  can  be  used  as  a  constraint  or  objective  function  in  a  parameter  optimization  problem.3 

The  process  illustrated  in  Figure  3  is  expressed  in  the  set  notation  expression: 

Clx  =  (^jRSEAnAS.  (1) 


7 


(j)  RSte  (k)  AS  (1)  Cix 

Figure  3:  Single  Coverage  Illustration  6  Satellite  Constellation 


4.3  Double  Coverage 

The  general  approach  to  the  double  coverage  case  is  similar  to  the  single  coverage  case  in  that  the  regions  of 
double  coverage  are  joined  by  union  operations  to  form  a  total  effective  coverage  region  that  is  intersected 
with  the  region  of  interest.  In  the  single  coverage  case,  each  effective  sensor  coverage  region  itself  exhibits 
the  desired  coverage  multiplicity.  In  contrast,  regions  of  double  coverage  are  determined  by  performing 
intersection  operations  between  unique  and  valid  pairs  of  effective  coverage  regions.  The  pairs  must  be 
unique  in  the  sense  that,  for  example,  comparing  satellites  A  and  B  for  overlap  is  equivalent  to  checking 
satellites  B  and  A  (i.e.  avoid  permutations  of  pairs  that  have  already  been  checked  because  they  only 
offer  redundant  information).  Similarly,  analyzing  for  overlap  between  satellites  A  and  A  is  also  an  invalid 
comparison  because  a  satellite  cannot  provide  double  coverage  on  its  own. 
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Figure  4  illustrates  the  sequence  of  Boolean  operations  necessary  to  determine  the  region  of  double 
coverage  in  the  region  of  interest,  C 2X-  Figures  4a-4c  and  Figures  4d-4f  illustrate  the  identification  by 
Boolean  intersection  of  two  unique  regions  exhibiting  double  coverage  between  neighboring  satellites.  These 
two  unique  regions  are  then  combined  using  a  union  operation,  as  shown  in  Figures  4g-4i.  As  more  regions 
of  double  coverage  are  identified,  they  are  also  joined  until  the  total  effective  ATH  coverage  region,  RSte, 
is  identified,  shown  in  Figure  4j.  RSte  is  then  intersected  with  AS  to  form  C2X,  the  region  of  double  ATH 
coverage  within  the  region  of  interest,  as  shown  in  Figure  41.  The  upper  bound  on  the  number  of  Boolean 
operations  involved  in  this  process  is  (n2  —  n)/ 2. 3 

The  sequence  of  Boolean  operations  can  be  expressed  in  set  notation  as: 

(n— 1  n  \ 

u  u  (RSEiinRSEi2)\  HAS,  (2) 

H=l*2=*l+1  / 

where  the  indices  on  the  finite  union  operators  are  specifically  chosen  to  avoid  invalid  or  redundant  satellite 
pair  comparisons.3 
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n 


(a)  RSe j 


(d)  RSE2 


(B)Ti 


(j)  RSte 


(b)  RSe2 


(k)  AS 


(c  )Ti 


(0  ^2 


(i)r3 


(1)  C-2X 


Figure  4:  Double  Coverage  Illustration  -  6  Satellite  Constellation 


4.4  Arbitrary  Coverage  Multiplicity 

Investigation  into  the  single,  double,  and  triple3  coverage  cases  leads  to  a  generalized  expression  to  efficiently 
identify  regions  of  any  desired  coverage  multiplicity,  denoted  by  p  (integer).  The  sequence  of  Boolean 
operations  is  compactly  described  in  set  notation  as: 

(n—p+l  n— p+2  n— 1  n  p  \ 

u  u  •  u  u  n^jcdis.  (3) 

ii  =  l  i2=il  +  l  ip— i=ip_ 2  +  1  ip=ip  — i  +  l  j=l  / 

The  primary  complexity  in  developing  this  expression  is  in  determining  the  appropriate  indices  on  the  finite 
union  operators  to  avoid  redundant  or  invalid  p-tuplets.3  For  values  of  p  =  1,2,  Equation  3  simplifies  to  the 
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cases  shown  in  Equations  1  and  2. 


4.5  Implementation 

The  set  notation  expressions  developed  in  Sections  4. 2-4.4  are  completely  implemented  in  both  MATLAB 
and  C-| — E  The  development  process  and  resulting  algorithms  are  extensively  documented,3  with  an  emphasis 
on  analysis  of  the  accuracy  of  discretized  region  approximation,  and  algorithm  performance.  Using  those 
results,  given  a  required  relative  accuracy  of  the  discretized  region,  an  appropriate  polygon  resolution  may  be 
selected  and  algorithm  runtime  can  be  estimated  for  both  the  MATLAB  and  CH — h  implementations  (runtime 
increases  linearly  with  polygon  resolution).3 

The  MATLAB  implementation  utilizes  the  General  Polygon  Clipping  Library,7  written  in  C.  Via  a 
MEX  gateway,  the  MATLAB  implementation  performs  polygon  clipping  operations  using  the  compiled 
binary.  Because  polygon  clipping  is  by  far  the  most  computationally  intensive  component  of  the  analysis, 
performance  of  the  MATLAB  implementation  is  only  slightly  slower  than  its  C++  counterpart.  The  C++ 
implementation  is  found  to  operate  approximately  4  times  faster  than  the  MATLAB  implementation,  and 
should  be  used  for  any  large  scale  problems  that  require  a  large  number  of  function  evaluations  to  solve.3 
However,  for  most  investigations,  that  are  typically  smaller  in  scale,  the  MATLAB  implementation  is  used 
exclusively  for  convenient  access  to  the  available  MATLAB  visualization  capabilities.  In  contrast  to  the 
MATLAB/C  implementation,  the  C++  implementation  has  yet  to  be  extended  to  the  volumetric  case,  as 
will  be  discussed  in  Section  6 
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5  Planar  Analysis  Example  Problems 

The  examples  presented  in  this  section  demonstrate  the  application  of  the  planar  constellation  analysis 
methodology  developed  in  the  previous  section  to  simple  constellation  optimization  problems.  Simple  time- 
invariant  models  are  optimized  first,  followed  by  several  time- varying  models  of  greater  complexity. 

5.1  Time-Invariant  Problems  —  Grid  Search 

Cases  where  the  distances  between  the  satellites  and  target  regions  remain  fixed  may  be  considered  time- 
invariant.  That  is,  the  amount  of  ATH  coverage  remains  constant  as  the  satellites  in  the  constellation  evolve 
along  their  orbit.  Based  on  this  assumption,  the  following  examples  demonstrate  the  use  of  the  discussed 
methodology  in  optimal  constellation  design.  The  first  example  presents  the  minimum  number  of  satellites 
required,  over  a  range  of  constellation  altitudes,  to  achieve  99.9%  single,  double,  and  triple  coverage  over 
a  range  of  altitudes.  The  second  example  expands  on  this  by  also  considering  the  omni-directional  sensor 
range  as  a  design  parameter. 

5.1.1  Example  1:  A  Single  Independent  Variable  Case 

A  planar  constellation  providing  ATH  coverage  to  an  Earth-centered  annular  target  region  is  considered.  The 
n  omni-directional  sensor  platforms  are  equally  distributed  in  a  single  circular  orbit.  Constellation  altitudes 
between  100  and  6000  km  are  considered  at  1  km  resolution.  The  remaining  problem  parameters  are  defined 
in  Table  1.  At  each  altitude,  the  minimum  constellation  populations  providing  at  least  99.9%  single,  double, 
or  triple  coverage  are  determined  with  a  simple  grid  search,  the  results  of  which  are  shown  in  Figure  5  and 
Table  2.  Just  as  coverage  area  is  computed  numerically,  so  is  the  area  of  the  target  region.  Due  to  roundoff 
and  truncation  error,  the  two  computed  areas  may  not  be  the  same,  despite  representing  the  exact  same 
regions.  To  prevent  this  from  causing  an  erroneous  result,  99.9%  coverage  is  considered  rather  than  100%. 
The  minimum-altitude  constellation  configurations  providing  single,  double,  and  triple  coverage  are  shown 
in  Figures  6a-6c. 


Constellation  Altitude  (km) 


Figure  5:  Minimum  Constellation  Population  vs.  Altitude 
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(a)  lx  cov.,  3  sat.,  h  =  999  km 


(b)  2x  cov.,  6  sat.,  h  =  1058  km 


(c)  3x  cov.,  10  sat.,  h  =  705  km 


Figure  6:  Smallest  Constellations  Providing  at  Least  99.9%  ATH  Coverage  at  Different  Multiplicities 


Table  1:  Example  1  Parameters 


Parameter 

Description 

Value 

ht 

tangent  height 

100  km 

hi 

lower  target  altitude 

1000  km 

hu 

upper  target  altitude 

5000  km 

R 

omni-directional  sensor  range 

10000  km 

m 

initial  polygon  resolution 

100  PPC 

Table  2:  Example  1  Optimal  Solutions 


Coverage  Mult. 

Optimal  Pop. 

Altitude  Range 

lx 

3  satellites 

999  -  1210  km 

2x 

6  satellites 

1058  -  1187  km 

3x 

10  satellites 

705  -  1113  km 
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5.1.2  Example  2:  A  Two  Independent  Variable  Case 


Expanding  upon  the  problem  in  Example  1,  in  addition  to  varying  circular  orbit  altitude  (19500-36000 
km),  variation  in  omni-directional  sensor  range  (17000-30000  km)  is  considered  in  a  400x400  grid.  Typically, 
sensor  range  is  a  fixed  quantity  depending  upon  available  hardware.  However,  such  an  analysis  may  be  useful 
during  a  trade  study  to  determine  the  minimum  sensor  performance  required  to  achieve  coverage  subject 
to  other  constraints.  The  problem  parameters  are  shown  in  Table  3.  The  minimum  number  of  satellites 
required  to  achieve  at  least  99.9%  single  coverage  across  the  phase  space  is  shown  in  Figure  7. 


Figure  7:  Minimum  Constellation  Size  For  99.9%  Single  Coverage  vs.  Altitude  and  Sensor  Range 


Table  3:  Example  2  Parameters 


Parameter 

Description 

Value 

ht 

tangent  height 

100  km 

hi 

lower  target  altitude 

20000  km 

hu 

upper  target  altitude 

36000  km 

m 

initial  polygon  resolution 

100  PPC 

5.2  Time-Invariant  Problems  —  Non-Linear  Programming 

To  illustrate  the  utility  of  the  numerical  ATH  coverage  models,  several  example  problems  are  solved  using  a 
Non-Linear  Programming  solver  (MATLAB’s  fmincon ).8  Prior  to  this,  a  simple  continuously  differentiable 
financial  model  describing  constellation  deployment  cost  is  developed  for  use  as  both  an  objective  function 
and  a  constraint  function.  Example  3  considers  minimization  of  deployment  cost  subject  to  a  constraint 
on  minimum  acceptable  ATH  coverage.  Example  4  illustrates  the  use  of  an  arbitrary  sensor  profile  in  a 
constellation  design  problem,  maximizing  single  ATH  coverage  subject  to  a  constraint  on  deployment  cost.3 

5.2.1  Simple  Financial  Model 

The  design  problems  described  in  this  section  involve  constraining  or  minimizing  financial  cost,  while  si¬ 
multaneously  maximizing  or  constraining  (respectively)  some  amount  of  coverage  at  a  specified  coverage 
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multiplicity.  Models  estimating  program  cost  for  a  satellite  or  constellation  are  used  throughout  govern¬ 
ment  and  industry  (see  Wertz).9  These  models  are  highly  complex  with  a  staggering  amount  of  parameters 
involved.  Additionally,  these  exhaustive  models  can  behave  discontinuously  in  many  regions,  i.e.  when  a 
payload  of  several  satellites  becomes  too  massive  for  one  launch  vehicle,  necessitating  another. 

The  objective  of  the  examples  in  this  section  is  not  to  illustrate  these  complex  cost  estimation  relations, 
but  to  illustrate  how  the  numerical  ATH  coverage  models  briefly  described  in  section  4  can  be  used  as  part 
of  a  constellation  design  process.  Consequently,  a  much  simpler  financial  model  is  desirable,  one  that  is 
continuous,  easily  described  to  the  reader,  and  easily  implemented  by  the  investigator. 

The  model  developed  here  focuses  on  two  parameters  as  cost  drivers  behind  constellation  deployment  - 
the  mass  of  each  spacecraft,  and  their  circular  orbit  altitude.  The  current  study  focuses  on  constellations 
composed  of  satellites  arranged  in  a  single  circular  orbit.  The  time-invariant  nature  of  this  set  of  problems 
allows  them  to  be  far  more  tractable  using  the  limited  computational  resources  available. 

Constellation  Deployment  Cost  According  to  Gordon,10  the  mass  of  a  sensor  or  antenna  can  be  crudely 
described  as  proportional  to  the  square  of  its  design  range.  Using  this  relationship,  spacecraft  mass  can  be 
approximated  by 

m  =  a  +  bR2.  (4) 

The  parameter  a  is  the  base  satellite  mass  (i.e.  if  no  sensor  were  installed  at  all),  and  b  represents  the 
mass/sensor-range  coefficient.  The  total  cost  of  constellation  deployment  is  then  described  by 

T  =  n(e  +  cm  +  dm(h  —  hTe  f)2).  (5) 

A  linear  relationship  between  the  number  of  satellites,  n,  and  deployment  cost  T  is  assumed.  The  parameter 
e  represents  a  base  cost  for  each  launch  system,  c  represents  the  additional  cost  per  unit  mass  to  achieve  the 
reference  circular  orbit  at  an  altitude  of  hle f.  The  cost  of  increasing  satellite  altitude  (h)  from  the  reference 
altitude  is  assumed  to  increase  quadratically,  thus  the  coefficient  d  describes  the  cost  per  unit  mass  vs.  the 
square  of  the  increase  in  altitude.  This  assumption  of  quadratic  behavior  with  respect  to  altitude  variation 
is  made  to  create  a  more  interesting  deployment  cost  model,  i.e.  not  just  quadratic  in  R,  but  also  in  h.  In 
all  analyses  presented  here,  hie f  is  considered  to  be  the  lower  bound  on  permissible  satellite  altitude. 

Mass  and  altitude  are  the  originally  intended  cost  drivers;  however,  by  Equation  4,  R  is  used  in  place  of 
mass  as  a  model  parameter. 

Parameter  Values  All  but  one  of  the  parameters  used  in  the  deployment  cost  model  remain  the  same  in 
both  Examples  1  and  2.  However  Example  2  considers  a  satellite  sensor  profile  of  arbitrary  shape,  where 
the  mass/sensor-range  coefficient  is  reduced  by  a  factor  of  four.  This  assumption  is  made  because,  as  will  be 
seen  in  Section  5.2.4,  the  arbitrary  region  encompassed  in  RSal b  is  much  smaller  than  in  an  omni-directional 
case  using  the  same  value  of  R.  The  resulting  financial  costs  computed  using  these  parameters  are  not 
necessarily  correlated  to  real  world  dollars,  and  are  only  intended  as  a  means  of  comparison  between  the 
solutions  and  problems  presented  here.  They  are,  however,  selected  with  the  aim  of  obtaining  results  on  the 
same  order-of-magnitude  as  a  more  exhaustive  analysis.9  These  parameters  are  summarized  in  Table  4. 

Table  4:  Financial  Model  Parameters  (see  Equations  4,  5) 


Parameter 

Value 

Description 

a 

500  kg 

base  satellite  mass 

boD 

6  x  10-5  kg/km2 

mass/sens-range  coefficient  (OD  sensors) 

^arb 

1.5  x  10-5  kg/km2 

mass/sens-range  coefficient  (arb.  sensors) 

C 

1  x  10-4  $M/kg 

cost-per-kg  at  reference  altitude 

d 

1  x  10"9  $M/kg-km2 

cost-per-kg  vs.  squared  altitude  increase 

e 

10  $M 

base  launch  vehicle  cost 

^•ref 

188  km 

satellite  reference  altitude 
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5.2.2  Non-Linear  Programming 

The  examples  presented  in  this  section  are  simple  constellation  design  problems,  and  their  solutions  are  easily 
obtained  by  parameter  optimization.  The  problems  are  parameterized  in  terms  of  a  finite  set  of  variables 
that  fully  define  a  unique  state  of  the  system.  The  solution,  or  optimal  values  of  these  variables  create  some 
optimal  condition  of  the  system  in  whatever  sense  the  investigator  defines  (by  finding  the  states  providing 
minimal  or  maximal  values  of  some  performance  index).  Although  the  examples  in  this  section  are  simple 
enough  to  be  analyzed  based  on  clever  plots  alone,  the  same  method  of  parameter  optimization  can  be  readily 
scaled  to  more  complex  problems,  for  which  visualization  of  the  phase  space  becomes  impossible. 

The  parameter  optimization  code  used  here  is  provided  in  the  MATLAB  Optimization  Toolbox  -  fmin- 
con ,8  which  is  a  package  of  several  different  non-linear  programming  (NLP)  packages.  Loosely  defined,  NLP 
is  a  process  by  which  some  objective  function  is  extremized  subject  to  equality  and/or  inequality  constraints, 
all  dependent  upon  a  finite  set  of  optimization  variables.  Additionally,  the  objective  function  and  the  con¬ 
straints  could  both  exhibit  non-linear  behavior  (NLP  can  be  considered  a  superset  of  linear  programming  - 
linear  programming  can  be  used  only  when  the  objective  and  constraint  functions  are  linear).  Specifically, 
the  sequential  quadratic  programming  (SQP)  algorithm  (see  Powell)11  in  fmincon  is  used  to  obtain  solutions 
to  the  example  problems  in  this  section.  The  general  problem  solved  by  NLP  is  as  follows: 

minimize  J  =  F(xp), 
subject  to  ceq(xp)  =  0, 
c(xp)  <  0, 

keeping  in  mind  that  maximizing  a  function  is  equivalent  to  minimizing  its  negative. 

5.2.3  Example  3  —  Min.  Budget  with  an  Area  Constraint 

The  constellation  designed  in  this  example  must  provide  single  coverage  to  at  least  99.9%  of  the  total 
area  enclosed  in  the  dual-altitude  band  shell.  This  value,  rather  than  100%,  is  arbitrarily  chosen  to  avoid 
numerical  issues  associated  with  roundoff  and  truncation  error  (however  small  they  may  be).  The  parameters 
describing  the  tangent  height  shell  (THS)  and  the  dual-altitude  band  shell,  AS/  are  shown  in  Table  5  along 
with  the  computed  area  inside  AS. 


(6) 

(7) 

(8) 


Table  5:  Example  3  Parameters 


Parameter 

Value 

Description 

Re 

6378.14  km 

assumed  Earth  radius 

ht 

100  km 

tangent  height  altitude 

h 

1000  km 

lower  target  altitude 

hu 

5000  km 

upper  target  altitude 

m 

100  PPC 

initial  polygon  resolution 

Aas 

235,307,769  W 

total  area  in  AS1 

Assume  that  the  satellites  are  equally  distributed  in  a  single  circular  orbit,  and  are  equipped  with  omni¬ 
directional  sensors.  Under  these  assumptions,  the  satellite  sensor  range,  R ,  and  the  satellite  altitude,  h,  are 
sufficient  to  fully  define  a  unique  state  of  the  system.  The  two-dimensional  set  of  variables  populates  the 
parameter  vector 

xp  =  [R  h]  .  (9) 

The  performance  index  that  is  minimized  with  respect  to  xp  (deployment  cost,  according  to  the  model 
developed  in  Section  5.2.1)  is  given  by 

J  =  T(R,h,n).  (10) 

Additionally,  three  inequality  constraints  are  present  as  follows: 


Amin  Aix 

-R 

/i-ref  h 


<  0. 


(11) 
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where  the  inequality  is  understood  to  apply  element-wise  between  vectors  c  and  0. 

The  parameters  and  solution  to  the  3  satellite  case  are  shown  in  Table  6. 

Figure  8a  shows  the  area  constraint  (the  solution  must  lie  inside  or  on  the  dotted  boundary).  From 
inspection  of  Figure  8c,  it  is  clear  to  see  that  the  converged  solution  does  correspond  to  a  constrained 
minimum  in  deployment  cost.  The  configuration  of  the  converged  solution,  shown  in  Figure  9b,  illustrates 
that  the  sensor  range  corresponds  to  the  intersection  between  neighboring  satellite  sensor  range  shells  and 
the  upper  target  altitude  shell  (UTAS),  as  expected  by  intuition. 

Table  6:  Example  3  -  Parameters  and  Solution 


Parameter 

Value 

Description 

n 

3 

number  of  satellites 

A 

-^min 

0.999Aas 

area  constraint 

Ro 

1000  km 

initial  guess,  R 

ho 

1000  km 

initial  guess,  h 

-*k)pt 

9889.940  km 

converged  cost-optimal  R 

^opt 

1038.909  km 

converged  cost-optimal  h 

ropt 

234.893  $M 

converged  cost  at  solution 

^lXopt  ^min 

-1.095  x  10-2  km2 

coverage  area  surplus 

5.2.4  Example  4  —  Arbitrary  Sensor  Profile 

This  problem  illustrates  how  arbitrary  sensor  profiles  can  be  used  with  the  numerical  ATH  coverage  models 
developed  in  Section  4  to  design  constellations.  The  satellite  sensor  cross-section  considered  here  is  by  no 
means  based  on  any  real-world  hardware.  Prior  to  addressing  the  actual  design  problem,  a  brief  discussion 
is  presented  that  illustrates  how  the  effective  satellite  range  shell,  RSe  is  obtained,  starting  from  a  simple 
sketch. 


Defining  the  Arbitrary  Effective  Range  Shell  One  of  the  fundamental  advantages  of  the  numerical 
analysis  of  ATH  coverage  is  that  investigations  are  not  restricted  simply  to  the  omni-directional  sensor 
case,  or  any  case  where  an  analytical  geometric  representation  of  the  sensor  cross-section  can  be  derived. 
Consider  the  case  shown  in  Figure  10.  This  sketch  represents  half  of  an  in-plane  sensor  cross-section,  and 
was  manually  drawn  by  the  investigator.  There  are  no  exact  mathematical  equations  governing  the  path, 
making  it  ‘arbitrary’  for  the  purposes  of  this  discussion. 

Upon  importation  of  the  drawing  via  digital  scanner,  it  is  slightly  adjusted  using  a  graphics  program 
to  smooth  the  boundary  and  correct  the  angle  at  the  end  of  the  curve,  making  it  perpendicular  to  the  x 
direction  (as  intended).  The  resulting  image  is  then  analyzed  using  a  figure  analysis  software  package,12 
producing  a  series  of  (x,y)  coordinates  describing  the  half-curve.  In  MATLAB,  the  half-curve  is  imported, 
the  sharp  tip  is  moved  to  the  origin,  and  the  lower  half  is  generated  by  replicating  the  ( x ,  y)  coordinates  in 
reverse,  while  changing  the  y  coordinates  to  their  negative  reciprocals.  The  coordinates  are  then  scaled  to 
produce  a  teardrop-shaped  sensor  cross-section  that  is  of  unit  length,  with  a  maximum  width  of  0.6.  The 
sharp  tip  is  centered  on  the  spacecraft,  as  shown  in  Figure  11.  This  curve  can  then  be  scaled  to  any  length, 
R ,  simply  by  multiplying  all  (x,y)  coordinates  by  R  (prior  to  any  coordinate  translation,  of  course). 

Further,  an  in-plane  attitude  of  the  spacecraft  (which  is  now  a  parameter,  as  the  sensor  is  no  longer 
omni-directional),  a  can  be  defined.  This  angle  is  defined  relative  to  the  plane  that  is  orthogonal  to  the 
satellite  position  vector.  Traditionally,  this  plane  is  referred  to  as  the  local  horizon.  However,  in  this  study, 
‘horizon’  refers  to  a  circular  boundary  around  the  Earth,  and  the  two  concepts  should  not  be  confused. 
In-plane  attitude  is  shown  in  Figure  12.  Due  to  the  symmetry  about  the  satellite  look-axis,  denoted  by  the 
arrow  on  the  satellite  body  in  Figure  12,  in-plane  attitude  is  constrained  by 


7 r  7T 

—  —  <  Q!  <  — . 
2  ~  ~  2 


(12) 


Once  scaled,  translated,  and  rotated  as  necessary,  this  arbitrarily  shaped  cross-section  represents  the  sensor 
range  shell,  RSarb- 
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(a)  Constraint  Function  (Single  Coverage  Area) 
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(b)  Constraint  Function  (Single  Coverage  Area) 
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Figure  8:  Example  3  -  Constraint  and  Objective  Contours 


18 


(a)  Initial  Guess 


(b)  Converged  Solution 

Figure  9:  Example  3  Initial  Guess  and  Converged  Solution 
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Figure  11:  Example  4  -  Arbitrary  Sensor  Cross-Section  Imported  to  MATLAB 
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Although  this  arbitrary  sensor  cross-section  represents  a  tear-drop  shape,  any  shape  could  be  used,  and 
placed  anywhere  in  relation  to  the  spacecraft.  The  only  caveat  is  that  the  sensor  region  must  be  assumed  to 
possess  some  form  of  symmetry  across  the  orbit-plane  such  that  ATH  coverage  measured  in-plane  somehow 
corresponds  to  ATH  coverage  in  a  three-dimensional  volumetric  sense. 

Having  defined  the  arbitrary  range  shell,  RSar b,  the  region  corresponding  to  the  tangent  height  triangle 
(THT)  is  removed  in  order  to  form  the  arbitrary  effective  range  shell,  RSe-  This  involves  the  formation  of 
a  polygon  representing  the  THT,  that  is  used  to  perform  the  polygon  clipping  operation  RS  —  THT  =  RSe- 
The  two  far-side  vertices  of  the  THT  can  be  determined  by  simple  geometry,  and  the  third  vertex  is  coincident 
with  the  satellite.  The  clipping  operation  is  illustrated  in  Figure  13,  where  the  THT  is  represented  by  the 
dashed  triangle,  and  the  shaded  region  represents  the  effective  range  shell,  RSe-  The  extent  of  the  THT 
need  only  exceed  the  maximum  extent  of  RSaib  to  ensure  there  are  no  leftover  regions  beyond  the  THT 
after  the  clipping  operation.  In  this  example,  the  distance  to  the  far  side  of  the  THT  (along  the  look-axis) 
is  chosen  to  be  1.1 1?. 

Possessing  the  means  to  generate  the  arbitrary  effective  range  shells,  it  is  then  straightforward  to  utilize 
the  numerical  ATH  coverage  models  developed  in  Section  4  to  evaluate  the  coverage  they  provide. 

Problem  Statement  and  Solution  The  problem  solved  in  this  example  is  maximization  of  single  ATH 
coverage  subject  to  a  deployment  constraint  for  an  arbitrarily  shaped  sensor  cross-section.  It  is  important  to 
note  that  the  financial  model  parameters  are  slightly  modified  for  this  example,  as  shown  in  Table  4,  where 
&arb  is  used.  This  parameter  is  adjusted  from  &od  to  reflect  the  reduced  coverage  area  relative  to  the  sensor 
range,  R. 

First,  the  parameter  vector  is  formed,  similarly  to  Example  1,  but  with  the  addition  of  the  in-plane 
attitude,  a  (uniform  across  the  constellation): 

xp  =  [R  h  a).  (13) 

Next,  the  objective  function  is  defined  as 

J  =  -Ai x  (i?,  h,  a,  n ),  (14) 

which  is  to  be  minimized  (to  maximize  single  coverage).  This  problem  is  subject  to  five  inequality  constraints 
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Figure  13:  Example  4  -  Clipping  the  THT  From  RSaio  to  Form  RSe 


that  are  written  in  vector  form  as 


c  = 


T(R,h,n)  —  Fmax 
-R 

href  h 
a  —  § 


<  0. 


(15) 


Note  that  the  financial  model,  described  in  Section  5.2.1,  does  not  depend  upon  in-plane  attitude,  a. 

The  parameters  used  in  this  problem  are  shown  in  Table  7.  The  initial  guess  is  intentionally  chosen  to 
be  far  away  in  all  three  variables  from  the  expected  (by  intuition)  solution.  Also  note  that  although  100 
PPC  is  the  polygon  resolution  used  for  the  altitude  shell,  the  arbitrary  range  shell  as  generated  in  MATLAB 
consists  of  110  vertices.  This  is  purely  an  artifact  of  the  curve  generation  process,  and  an  order-of-magnitude 
similarity  is  considered  acceptable. 


Table  7:  Example  4  -  Parameters 


Parameter 

Value 

Description 

Re 

6378.14  km 

assumed  Earth  radius 

ht 

100  km 

tangent  height  altitude 

hi 

1000  km 

lower  target  altitude 

fa'll 

5000  km 

upper  target  altitude 

m 

100  PPC 

initial  polygon  resolution 

n 

7 

number  of  satellites 

r 

1  max 

250  $M 

deployment  cost  constraint 

Ro 

1000  km 

initial  guess,  R 

ho 

4000  km 

initial  guess,  h 

01 0 

45° 

initial  guess,  a 
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After  41  SQP  iterations,  the  optimizer  converges  to  the  solution  shown  in  Table  8.  In  order  to  visualize 
this  three  degree-of-freedom  problem,  contours  at  the  solution  are  produced  holding  one  variable  fixed, 
creating  a  pair  of  contour  plots  for  ATH  coverage  area  and  deployment  cost.  Figure  14  depicts  overall 
and  detail  contours  when  holding  in-plane  attitude  constant  at  the  solution.  Note  that  the  detail  contour 
(Figure  14b)  uses  a  rescaled  color  map,  and  the  colors  do  not  correspond  to  the  colors  in  the  overall  contour 
(Figure  14a).  These  figures  illustrate  that  coverage  area  only  decreases  in  any  feasible  direction  in  R  or 
h  (i.e.  not  across  the  constraint  boundary).  Similarly,  holding  h  fixed  at  the  solution  yields  the  contours 
shown  in  Figure  15.  Just  as  in  the  case  of  constant  in-plane  attitude,  any  variation  in  R  and  a  in  a  feasible 
direction  results  in  a  decrease  in  ATH  coverage  area.  As  with  earlier  examples  constrained  by  deployment 
cost,  because  of  the  quadratic  behavior  of  the  cost  function  in  both  R  and  h,  the  solution  must  fall  on  or  to 
the  left  of  the  constraint  boundary. 

While  these  contour  illustrations  are  by  no  means  sufficient  condition  to  declare  optimality  of  the  solution, 
the  SQP  optimizer  in  fmincon  does  return  a  flag  indicating  that  it  has  found  a  constrained  locally  optimal 
solution.  The  MATLAB  SQP  implementation  utilizes  a  variation  of  the  Karush-Kuhn- Tucker  (KKT)  condi¬ 
tions13  as  one  of  several  criteria  to  determine  optimality.  Generally,  KKT  is  only  a  necessary  condition  for 
optimality.  However,  for  continuously  differentiable  and  convex  objective  and  constraint  functions  (as  is  the 
case  in  this  problem),  KKT  also  provides  sufficient  condition  for  optimality.14 

Table  8:  Example  4  -  Solution 


Parameter 

Value 

Description 

^opt 

11396.402  km 

converged  area-optimal  R 

^opt 

897.559  km 

converged  area-optimal  h 

C^opt 

-5.937° 

converged  area-optimal  a 

Xopt 

214,239,724  km2 

single  coverage  area  at  solution 

r  opt  r  max 

0  $M 

budget  overrun 
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5.3  Time-Varying  Problems 

In  cases  with  elliptical  orbits,  more  complex  target  regions,  or  even  constellations  with  populations  distributed 
across  multiple  circular  and/or  elliptical  orbits  (in  the  same  plane),  the  ATH  coverage  amount  is  time- varying 
in  general.  The  overlap  between  sensor  coverage  and  target  regions  varies  continuously. 

Time-varying  problems  may  be  analyzed  using  the  techniques  developed  in  this  study  by  performing 
instantaneous  coverage  evaluations  at  specified  times  throughout  a  time  interval  of  interest.  Figure  18  shows 
the  ATH  coverage  provided  by  a  planar  constellation  of  6  sensor  platforms  in  arbitrarily  prescribed  orbits 
with  arbitrarily  defined  sensor  profiles  over  a  prescribed  time  interval.  A  sufficiently  long  time  interval,  and 
a  sufficiently  short  time-step  must  be  selected  to  ensure  the  analysis  adequately  characterizes  the  behavior 
of  the  system.  In  a  system  with  periodic  behavior,  an  appropriate  time  interval  is  one  period.  However, 
non-periodic  cases  such  as  that  shown  in  Figure  18  are  more  difficult  to  analyze  because  the  judgment  of 
the  analyst  must  dictate  an  appropriate  time  interval.  The  case  shown  in  Figure  18,  and  animated  in  the 
accompanying  video  file,  features  6  satellites  with  arbitrary  sensor  profiles  in  arbitrarily  chosen  coplanar 
elliptical  orbits.  Each  sensor  platform  maintains  a  fixed  (and  prescribed)  attitude  with  respect  to  the  local 
horizon  at  all  times.  The  target  region  is  an  Earth-fixed  dual-altitude  band  region  bounded  in  longitude. 

Using  this  numerical  approach,  constellation  design  problems  can  be  addressed  using  various  parameter 
optimization  techniques. 

5.3.1  Example  5  Analysis  of  Coverage  by  Elliptical  Orbits 

As  in  the  time-invariant  problems,  an  annular  target  region  is  considered.  However,  due  to  the  time- 
invariant  nature  of  the  target  region  itself,  a  time-invariant  solution  is  optimal,  as  is  demonstrated  in  this 
example. 

The  problem  is  posed  as  a  Mixed-Integer  Non-Linear  Programming  (MINLP)  problem  and  solved  using 
MIDACO.15  MIDACO  is  a  zeroth  order  heuristic  solver  and  uses  an  ant  colony  optimization  (ACO)  ap¬ 
proach.16  Because  ACO  is  a  heuristic  approach  to  MINLP,  there  are  no  analytical  optimality  criteria  for 
non-convex  problems  such  as  this,  thus  the  algorithm  is  allowed  to  run  until  it  ceases  improvement  upon  the 
solution. 

A  constellation  composed  of  two  elliptical  orbits  with  opposite  periapsis  directions  is  considered.  Each 
orbit  is  initially  populated  by  three  satellites,  each  group  with  dissimilar  omni-directional  sensor  performance. 
Within  each  orbit  the  satellites  are  equally  spaced  in  mean  anomaly.  The  objective  is  the  minimize  the 
total  population  of  the  constellation,  while  ensuring  continuous  single  coverage  of  the  target  region.  The 
formulation  allows  for  one  or  both  orbits  to  have  zero  population  (although  zero  population  in  both  orbits 
results  in  an  obvious  violation  of  the  continuous  coverage  constraint).  Problem  parameters  are  summarized 
in  Table  9. 


Table  9:  Example  5  -  Parameters 


Parameter 

Description 

Value 

ht 

tangent  height 

100  km 

hi 

lower  target  altitude 

1000  km 

hu 

upper  target  altitude 

10000  km 

k 

number  of  distinct  orbits  in  constellation 

2  orbits 

Ri 

omni-directional  sensor  range  of  sats.  in  orbit  1 

5000  km 

R2 

omni-directional  sensor  range  of  sats.  in  orbit  2 

10000  km 

m 

initial  polygon  resolution 

100  PPC 

Table  10  shows  the  start  point  and  solution  after  300  seconds  (682  function  evaluations).  From  an 
infeasible  start  point  (coverage  gaps),  as  seen  in  Figure  19-a,  a  feasible  12  satellite  state  is  identified  by 
the  15th  function  evaluation.  By  the  98th  function  evaluation,  the  solution  presented  in  Figure  19-b  is 
converged  upon.  With  6  satellites  distributed  in  a  single  circular  orbit,  the  solution  is  time- invariant  when 
only  concerned  with  the  quantity  of  coverage,  as  expected. 
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Figure  14:  Example  4  -  Objective  Function  (Single  Coverage  Area,  Constant  a  =  —5.937°  at  Solution) 
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Figure  15:  Example  4  -  Objective  Function  (Single  Coverage  Area,  Constant  h  =  897.559  km  at  Solution) 
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Figure  16:  Example  4  -  Constraint  Function  (Deployment  Cost) 
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(a)  Initial  Guess 
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Figure  18:  ATH  coverage  over  time  -  arbitrary  sensor  profiles,  asymmetric  target  region  (see  accompanying 
video  #1) 


Table  10:  Example  5  -  Start  Point  &  Solution  (682  func.  evals,  300s) 


Decision  Variable 

Description 

Start  Point 

Solution 

a 

semi-major  axis  (both  orbits) 

10000  km 

8588.0 

e 

eccentricity  (both  orbits) 

0.25 

0 

M20 

M  at  epoch  for  lead  sat.  in  orbit  2 

0  rad 

0  rad 

ni 

orbit  1  population 

3  satellites 

0  satellites 

n2 

orbit  2  population 

3  satellites 

6  satellites 
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Figure  19:  Example  5  -  Start  point  to  solution  with  continuous  lx  coverage  constraint  (see  accompanying 
video  #2) 
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5.3.2  Example  6  —  Continuous  Observation  of  CONUS  GEO  Satellites 


In  contrast  to  the  previous  example,  an  asymmetric  Earth-fixed  target  region  is  considered  here.  From  a 
geostationary  orbit,  the  target  region  is  defined  as  ±1000  km  in  altitude,  between  148°W  and  61°W  longitude 
(populated  by  satellites  serving  the  continental  United  States). 

One  immediately  obvious  solution  corresponds  to  placing  the  satellites  on  a  geostationary  orbit.  However, 
because  this  particular  arrangement  is  itself  time- invariant,  it  is  not  particularly  useful  in  demonstrating 
the  algorithm’s  capability  for  addressing  time- varying  problems.  Instead,  consider  an  alternate  arrangement 
where  the  constellation  is  composed  of  satellites  placed  across  four  smaller  identical  orbits,  equally  distributed 
in  periapsis  direction.  Unlike  previous  examples,  the  satellites  are  no  longer  equally  distributed  within  each 
orbit,  and  are  instead  equally  distributed  within  a  range  in  mean  anomaly.  The  positioning  (mean  anomaly 
at  epoch)  of  each  group  is  prescribed  such  that  apoapsis  of  the  center  of  each  satellite  group  occurs  as  the 
target  region  is  centered  above  apoapsis  of  each  orbit.  In  order  to  maintain  this  synchronization,  the  time 
period  (and  thus  semi-major  axis)  of  the  orbits  is  prescribed  to  revisit  geostationary  altitude  an  integer 
number  of  times  per  day  (twice  daily  in  this  example). 

The  objective  in  this  example  is  to  minimize  the  satellite  sensor  range,  i.e.  determine  the  smallest  sensor 
range  capable  of  covering  the  target  region  continuously.  The  problem  is  subsequently  posed  as  a  non¬ 
linear  programming  (NLP)  problem,  where  all  integer  parameters  are  prescribed,  and  approached  using  the 
interior-point  solver  in  fmincon .8  The  decision  variables  are  satellite  sensor  range,  R ,  orbit  eccentricity  e, 
and  satellite  group  spread,  AM.  Problem  parameters  are  summarized  in  Table  11. 

Table  11:  Example  6  Parameters 


Parameter 

Description 

Value 

ht 

tangent  height 

100  km 

hi 

lower  target  altitude 

41164.13  km 

hu 

upper  target  altitude 

43164.13  km 

A  i 

western  target  longitude 

148°W 

^ U 

eastern  target  longitude 

61°W 

k 

number  of  distinct  orbits  in  constellation 

4  orbits 

n 

total  constellation  population 

8  satellites 

m 

initial  polygon  resolution 

100  PPC 

Table  12:  Example  6  -  Start  Point  &  Solution 


Decision  Variable 

Description 

Start  Point 

Solution 

R 

omni-directional  sensor  range 

12000  km 

17444  km 

e 

orbit  eccentricity 

0.50 

0.40 

AM 

group  spread 

90°  km 

89.4° 

Table  12  shows  the  infeasible  start  point  and  solution,  obtained  after  32  iterations.  The  start  point  and 
solution  constellations  are  shown  in  Figures  20a  and  20b  respectively.  Although  the  number  of  satellites  in 
each  orbit  and  the  constellation  as  a  whole  remain  fixed,  this  example  demonstrates  how  the  numerical  ATH 
coverage  model  may  be  used  in  an  NLP-driven  constellation  design  process. 
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(b)  Solution 


(a)  Start  Point 
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Figure  20:  Example  6  -  Start  point  to  solution  with  continuous  1  x  coverage  constraint  (see  accompanying 
video  #2) 
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6  Non-Planar  Constellation  Analysis 

Volumetric  coverage  analysis  for  potentially  non-planar  constellations  is  accomplished  by  extending  the 
previously  described  planar  analysis  code.  When  considering  satellites  in  non-coplanar  orbits,  there  is  no 
single  cross-sectional  plane  to  compute  planar  coverage  area  within.  Consequently,  the  coverage  area  in  any 
one  plane  has  no  correlation  to  the  coverage  volume  of  a  given  non-coplanar  constellation.  To  directly  analyze 
volumetric  coverage,  the  planar  analysis  code  analyzes  coverage  area  of  sensor  and  target  region  cross-sections 
within  a  finite  set  of  equally  spaced  planes  offset  from  a  prescribed  analysis  plane  (i.e.  the  xy  plane  in  an 
Earth-centered  inertial  reference  frame) .  When  coverage  area  in  each  plane  is  determined,  it  is  multiplied  by 
the  distance  between  the  offset  planes,  resulting  in  an  approximate  volume  for  a  particular  ‘cutting  plane’  of 
constellation  coverage.  The  sum  of  these  ‘cutting  plane’  volumes  is  the  approximate  coverage  volume  for  the 
entire  constellation.  The  accuracy  of  the  coverage  volume  approximation  is  directly  proportional  to  ‘cutting 
plane’  density  over  the  region  covered  by  the  constellation. 

The  fundamental  challenge  in  analyzing  volumetric  coverage  of  constellations  then  becomes  cross-section 
determination  of  sensor  and  target  regions.  Omni-directional  sensor  profiles  are  considered  in  the  present 
investigation,  although  this  approach  may  be  used  for  any  sensor  profile,  provided  a  means  exists  to  generate 
cross-sections  at  arbitrary  orientations  to  the  profile. 

In  this  section  the  methods  by  which  target  region  cross-sections  and  ATH  sensor  cross-sections  are 
obtained  for  the  omni-directional  sensor  case  are  developed.  These  results  are  applied  to  example  problems 
in  the  next  section,  demonstrating  volumetric  coverage  analysis  capability  for  the  omni-directional  sensor, 
non-coplanar  orbit  case. 
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6.1  ATH  Sensor  Coverage  Cross-Section  Generation 

Figure  21  shows  several  representative  cross-sections  of  an  approximated  omni-directional  sensor  coverage 
volume.  The  red  dot  indicates  the  location  of  the  satellite,  and  the  missing  section  is  the  tangent  height 
cone  (THC)  -  the  region  of  the  omni-directional  volume  that  is  below  the  horizon,  or  against  an  Earth 
background. 


(a)  Degenerate  Conic  (b)  Hyperbolic  Edge 


(c)  Disk  (d)  Overall 


Figure  21:  Offset  cross-sections  of  the  sensor  coverage  volume 


Figure  21a  shows  the  familiar  sensor  region  cross-section  seen  in  prior  planar  constellation  investigations. 
The  ‘cutting  plane’  intersects  the  position  of  the  satellite  and  the  axis  of  the  THC  (the  satellite  position 
vector),  creating  a  boundary  composed  of  a  circular  arc  and  a  degenerate  conic  section.  Figure  21b  shows 
a  cross-section  composed  of  a  circular  arc  and  a  hyperbolic  conic  section,  while  Figure  21c  shows  a  disk 
cross-section,  not  intersecting  the  THC. 

To  analyze  constellations  composed  of  satellites  in  non-coplanar  orbits,  sensor  coverage  volume  cross- 
sections  must  be  available  for  any  ‘cutting  plane’  angle,  (j).  The  cutting  plane  angle  refers  to  the  acute  angle 
between  the  satellite  position  vector  and  the  cross-section  plane  and  is  bounded  by  —  7t/2  <  (f>  <  7r/2.  Note 
that  4>  =  0  in  the  case  shown  in  Figure  21.  Figure  22  shows  a  more  general  case  with  a  positive  (f>,  featuring 
cross-sections  composed  of  circular  arcs  and  hyperbolic  conic  sections. 


(a)  Hyperbolic  Edge  (b)  Degenerate  Conic 


(c)  Hyperbolic  Edge 


(d)  Disk 


Figure  22:  Offset  cross-sections  of  the  sensor  coverage  volume  for  a  non-zero  cutting  plane  angle 


For  the  omni-directional  sensor  coverage  case,  given  an  omni-directional  sensor  range,  satellite  position, 
and  tangent  height,  vertices  outlining  cross-sections  of  the  ATH  sensor  coverage  volume  may  be  determined 
in  closed  form  by  solving  various  geometrically  obtained  equations. 

The  eccentricity  of  the  THC  conic  section  at  cutting  plane  (j)  may  be  computed  by  the  relation 


cos  \(j)\ 

COS"f 


(16) 


where  7  is  the  interior  half- angle  of  the  THC.  The  radius  of  the  base  of  the  THC  truncated  at  the  intersection 
of  the  THC  and  sensor  radius  (i.e.  the  ‘edge’  of  the  missing  region  in  Figures  21d  and  22e)  is  given  by 


to  =  R  sin  7. 


(17) 
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The  eccentricity  of  the  THC  conic  section  is  used  to  identify  different  cases  by  which  cross-sections  of  the 
omni-directional  ATH  coverage  region  are  determined.  For  0  <  e  <  1  —  ne,  the  cross-sections  are  either  disks, 
or  compound  regions  composed  of  circular  and  elliptical  arcs.  For  1  —  ne  <  e  <  1  +  ne,  the  cross-sections  are 
either  disks,  or  compound  regions  composed  of  circular  and  parabolic  arcs.  For  e  >  1  -I -ne,  the  cross-sections 
are  either  disks  or  compound  regions  composed  of  circular  and  hyperbolic  arcs.  Degenerate  conics  result 
when  the  cutting  plane  passes  through  the  satellite  position  (apex  of  the  THC).  This  occurs  when  the  cutting 
plane  offset,  h,  is  zero.  Note  that  in  the  implementation  degenerate  and  non-existent  cases  are  tested  for 
first  in  order  to  simplify  the  logical  structure  of  the  proceeding  code.  Thus,  the  degenerate  and  non-existent 
cases,  D  and  E  are  presented  prior  to  cases  A-C. 

The  conditions  for  a  parabolic  conic  section  of  the  THC  are  modified  with  respect  to  classical  analytical 
consideration  for  numerical  reasons.  In  practice,  on  a  finite  precision  computing  platform,  the  near-parabolic 
hyperbolic  or  elliptical  conic  sections  are  difficult  to  compute  accurately  due  to  the  discontinuity  in  several 
of  the  conic  section  parameters  in  the  neighborhood  of  e  =  1.  For  this  reason,  an  expanded  parabolic 
interval  of  ±ne  is  implemented,  where  near-parabolic  conic  sections  are  assumed  to  be  parabolic  and  are 
determined  by  relationships  derived  for  the  parabolic  case.  The  parameter  e  is  the  ‘machine  epsilon’  for  the 
computer  platform  being  used,  and  represents  the  distance  to  the  nearest  number  to  1  that  is  differentiable 
by  the  system  from  1.  n  is  a  parameter  that  is  determined  by  experimentation  to  encompass  the  region  of 
inaccuracy  surrounding  e  =  1,  and  is  assumed  to  be  1  x  108  in  the  current  study. 

The  rest  of  this  section  documents  the  geometric  relationships  necessary  to  define  the  omni-directional 
sensor  cross-sections  at  arbitrary  cutting  plane  angles.  This  capability  is  necessary  for  generalized  volumetric 
ATH  coverage  computation  using  the  volumetric  layering  approximation  method  investigated  in  this  study. 
The  cross-sections  are  derived  in  a  Cartesian  coordinate  system  ( xyz )  such  that  the  analysis  plane  is  parallel 
to  the  x  —  y  plane  in  the  general  problem,  the  z  axis  is  parallel  to  the  z  axis  in  the  general  problem,  and  the 
x  axis  points  toward  the  Earth  parallel  to  the  projection  of  the  satellite  position  vector  into  the  x  —  y  plane. 
Transformation  from  the  xyz  frame  back  to  the  xyz  frame  of  the  general  problem  is  discussed  at  the  end  of 
this  section. 

6.1.1  Case  D:  h  =  0 

When  the  cutting  plane  intersects  the  satellite  location  (the  apex  of  the  THC),  a  degenerate  conic  results. 
This  conic  section  is  either  a  single  point,  a  line  that  encompass  zero  area  (in  which  case  the  sensor  cross- 
section  is  regarded  as  a  disk) ,  or  a  triangular  region  bounded  by  two  lines  intersecting  at  the  apex  (resulting 
in  a  sensor  cross-section  resembling  the  in-plane  planar  case  cross-sections  of  previous  discussion).  Cutting 
plane  positions  leading  to  this  case  are  shown  in  Figure  23.  The  associated  cutting-plane  cross-sections  of 
the  ATH  coverage  volumes  in  Figure  23  are  shown  in  Figure  24. 


(a)  Dl:  7  >  \cj,\  (b)  D2:  7  =  \cj>\  (c)  D3:  7  <  |0| 


Figure  23:  Case  D 
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(a)  Dl:  7  >  |0| 


(b)  D2:  7  =  |0| 


(c)  D3:  7  <  |0| 


Figure  24:  Case  D 


Case  Dl:  7  >  \(j>\  The  degenerate  conic  section  is  bounded  by  two  lines  intersecting  at  the  satellite, 
separated  by  the  interior  half  angle 


7'  =  cos  1 


cos  7 
cos \<j>\ 


(18) 


The  circular  portion,  spanned  clockwise  is  bounded  by  the  angles 


%l>i  =  27T  -  7',  ijjf  =  7' 


(19) 


as  measured  in  the  cutting  plane  from  the  x  direction. 

Thus,  the  cross-section  is  a  disk  of  radius  R  missing  a  27'  wide  sector,  as  shown  in  Figure  24a. 


Cases  D2,  D3:  7  =  \(j>\  or  7  <  \cf>\  The  cutting  plane  does  not  intersect  the  interior  of  the  THC,  thus  the 
conic  section  of  the  THC  is  a  point  or  a  line  with  zero  area.  The  sensor  cross-section  is  assumed  with  no 
loss  of  accuracy  to  be  a  disk  of  radius  R ,  as  shown  in  Figures  24b  and  24c. 


6.1.2  Case  E:  \h\  >  R 

The  cutting  plane  lies  outside  the  radius  of  the  omni-directional  sensor  volume,  thus  the  cross-section  does 
not  exist.  Figure  25  shows  cutting  plane  offsets  (h)  leading  to  this  case.  The  associated  cross-sections  do 
not  exist,  thus,  they  are  not  shown. 


6.1.3  Case  Al:  e  >  1  +  ne  and  <f>>  0 

The  THC  cross-section  is  strictly  hyperbolic,  thus,  sensor  region  cross-sections  may  be  disks,  non-existent, 
or  defined  by  a  compound  boundary  between  a  circular  arc  and  a  hyperbola.  Excluding  the  cases  where  the 
cross-section  is  determined  to  be  a  disk,  the  problem  is  determining  the  hyperbolic  parameters  a  and  b  in 
the  canonical  form  of  a  hyperbola,  and  determining  the  center-point  of  the  conic-section,  s' .  s'  is  centered 
between  the  ‘first  point,’  pi ,  of  the  THC  conic  section,  and  the  ‘first  point’  of  the  imaginary  rear-half  of  the 
THC,  q\.  ‘First  point’  refers  to  the  closest  point  on  each  hyperbola  to  the  center.  By  finding  the  location 
of  a  ‘second  point,’  P2+,  where  the  conic  section  intersects  the  omni-directional  sensor  radius  in  the  cutting 
plane,  the  appropriate  angular  intervals  in  ip  may  be  determined  to  define  the  hyperbolic  and  circular  curves. 

For  the  hyperbolic  segments,  the  radius  of  the  hyperbola  with  respect  to  the  focus,  /,  located  on  the  x 
axis  at  xs'  +  ae,  is  given  by 

a(e2  —  1) 

^hyperbola  =  Z  7-  (20) 

l-ecosV’ 
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(a)  h  >  R 


(b)  h  <  -R 


Figure  25:  Case  E 


The  radius  of  the  outer  circular  arc,  bounding  the  omni-directional  sensor  radius  in  the  ATH  regions  is 
given  with  respect  to  the  origin  by 

^circular  arc  ~  \/ A7  .  (-  1  ) 

Case  A1A:  R  >  h  and  h  >  R  sin  (7  —  |0|)  The  cutting  plane  lies  above  the  upper  intersection  of  the  THC 
with  the  omni-directional  sensor  radius,  leading  to  a  disk-shaped  sensor  cross-section  with  radius  rc irCuiar  arc- 
Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  26.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  26  are  shown  in  Figure  27. 

Case  A1B:  I?  sin  (7  —  |^>|)  >  h  and  h  >  0  The  cutting  plane  intersects  the  THC  within  the  omni¬ 
directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between  a 
circular  arc  and  a  hyperbola. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  28.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  28  are  shown  in  Figure  29. 

First,  the  coordinates  of  the  ‘first  point’  on  the  THC,  p±,  are  determined  by 

_/  0  :7-|^l  =  f 

*”“t  sd™  f  ’  (22) 

Vpi  =  0- 

Similarly,  the  ‘first  point’  on  the  THC,  q±  is  determined  by 

_/  0  :  7+ |^|  =  | 

^  1  ^  +  1^?  ’  (23) 

Vqi  =  fl¬ 
it  can  be  shown  by  trigonometry  then  that  using  the  intermediate  values  q  and  g ,  found  by  the  relations 

,24, 
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z  z 


Figure  27:  Case  AlA 
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9=  \  R-  — 


h 


cos  7 

sm7  —  \<p\  J  cos\<j)\ 
the  coordinates  of  the  ‘second  point,’  P2+ ,  can  be  determined  by 

Xp2+  =  xPl  +  g , 


yP2+  =  \/m2  -  q2. 

It  then  follows  that  the  hyperbolic  parameters  a,  b,  and  xs>  may  be  determined  by 

Xr> i  Xn 


(25) 


(26) 


1 


b  =  a\J  e2  —  1, 


(27) 


(28) 


In  the  implementation,  the  hyperbolic  span  is  split  into  two  regions,  such  that  the  initial  point  of  the 
outer  polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of 
visualization  functions.  The  upper  half  of  the  hyperbola  is  thus  spanned  by  the  angles 


4>i  =  7T, 

ipf  =  arctan 


Up2+ 


XP2+  ~  (x's  +  ae)  J  ’ 


(29) 


relative  to  the  focus,  /.  The  radius  relative  to  /  is  determined  using  Equation  20.  The  circular  arc,  relative 
to  the  origin  is  spanned  by 


ipi  =  arctan 


(  VP2+ 


\x. 


P  2+ 


tpf  =  2tt  —  arctan 


(  yp2+_ 


V- 


'P  2  + 


(30) 


with  a  radius  given  by  Equation  21.  Finally,  the  lower  half  of  the  hyperbola  is  spanned  by 

VP2  + 


i pi  =  —  arctan 

lj)f  =  —7 T. 


bP  2+ 


-  (x’s  +ae)J  ’ 


(31) 


Case  A1C:  0  >  h  and  h  >  — JRsin (7  +  |0|)  The  cutting  plane  intersects  the  THC  within  the  omni¬ 
directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between  a 
circular  arc  and  a  hyperbola. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  30.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  30  are  shown  in  Figure  31. 

First,  the  coordinates  of  the  ‘first  point’  on  the  THC,  pi,  are  determined  by 


i-pi 


:7  +  H  =  | 
tan  (7+101)  :  7  +  \4>\  7^  f 


0 

1  h\ 


(32) 


Vpi  =  0. 

Similarly,  the  ‘first  point’  on  the  THC,  q\  is  determined  by 


Xni  = 


Vgi  =  °- 


:  7  —  101  =  f 

tan  (7— 10|)  :  7  -  I'/’l  0  f 


0 


(33) 
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It  can  be  shown  by  trigonometry  then  that  using  the  intermediate  values  g  and  q ,  found  by  the  relations 


9  ={R-  \Jh2  +  xh) 


9  = 


m  —  g- 


sm  ( 


cos  7 
cos\4>\ ' 

+  7) 


cos  7 


the  coordinates  of  the  ‘second  point,’  P2+,  can  be  determined  by 


2+ 


=  xPl  +  g, 


yP2+  =  V™-2  ~  Q2- 

It  then  follows  that  the  hyperbolic  parameters  a,  b,  and  xs'  may  be  determined  by 

Xr )i  Xa, 


(34) 

(35) 

(36) 


a  = 


b  =  a\J  e2  —  1, 


xs'  =  xv,  -  a. 


(37) 


(38) 


In  the  implementation,  the  hyperbolic  span  is  split  into  two  regions,  such  that  the  initial  point  of  the 
outer  polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of 
visualization  functions.  The  upper  half  of  the  hyperbola  is  thus  spanned  by  the  angles 


4>i  =  7T, 

ipf  =  arctan 


Up2+ 


Xp2+  -  (x's  +  ae)  J  ’ 


(39) 


relative  to  the  focus,  /.  The  radius  relative  to  /  is  determined  using  Equation  20.  The  circular  arc,  relative 
to  the  origin  is  spanned  by 


ipi  =  arctan 


f  VP2+ 


'V  2+ 


ipf  =  2-7T  —  arctan 


f  VP2+ 

V  XP2+ 


(40) 


with  a  radius  given  by  Equation  21.  Finally,  the  lower  half  of  the  hyperbola  is  spanned  by 

Up2+ 


ipi  =  —  arctan 

=  —  7T. 


bP  2  + 


-  (x'3  +  ae)  J  ’ 


(41) 


Case  AID:  —  A2  sin  ( | <^> |  +7 )  >  h  and  h  >  —R  The  cutting  plane  lies  below  the  lower  intersection  of  the 
THC  with  the  omni-directional  sensor  radius  leading  to  either  a  disk-shaped  sensor  cross-section  with  radius 
^circular  arc,  or  no  region  at  all  if  \<j>\  +  7  >  7r/2. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  32.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  32  are  shown  in  Figure  33. 

6.1.4  Case  A2:  e  >  1  +  ne  and  (f>  <  0 

The  THC  cross-section  is  strictly  hyperbolic,  thus,  sensor  region  cross-sections  may  be  disks,  non-existent, 
or  defined  by  a  compound  boundary  between  a  circular  arc  and  a  hyperbola.  Excluding  the  cases  where  the 
cross-section  is  determined  to  be  a  disk,  the  problem  is  determining  the  hyperbolic  parameters  a  and  b  in 
the  canonical  form  of  a  hyperbola,  and  determining  the  center-point  of  the  conic-section,  s' .  s'  is  centered 
between  the  ‘first  point,’  pi ,  of  the  THC  conic  section,  and  the  ‘first  point’  of  the  imaginary  rear-half  of  the 
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Figure  30:  Case  AlC 
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Figure  31:  Case  AlC 
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Figure  32:  Case  AID 
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Figure  33:  Case  AID 


43 


THC,  q\.  ‘First  point’  refers  to  the  closest  point  on  each  hyperbola  to  the  center.  By  finding  the  location 
of  a  ‘second  point,’  P2+ ,  where  the  conic  section  intersects  the  omni-directional  sensor  radius  in  the  cutting 
plane,  the  appropriate  angular  intervals  in  0  may  be  determined  to  define  the  hyperbolic  and  circular  curves. 

For  the  hyperbolic  segments,  the  radius  of  the  hyperbola  with  respect  to  the  focus,  /,  located  on  the  x 
axis  at  xs'  +  ae,  is  given  by 

°(e2  —  1)  ,  , 
^hyperbola  =  “  7  (42 ) 

1  —  e  cos  ip 

The  radius  of  the  outer  circular  arc,  bounding  the  omni-directional  sensor  radius  in  the  ATH  regions  is 
given  with  respect  to  the  origin  by 

^circular  arc  =  \J  R2  ~  h2 .  (43) 

Case  A2A:  R  >  h  and  h  >  i?sin(7  +  |</>|)  The  cutting  plane  lies  above  the  upper  intersection  of  the 
THC  with  the  omni-directional  sensor  radius  leading  to  either  a  disk-shaped  sensor  cross-section  with  radius 
^circular  arc,  or  no  region  at  all  if  101  +  7  >  71-/2. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  34.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  34  are  shown  in  Figure  35. 


Case  A2B:  f?sin  (7  +  |0|)  >  h  and  h  >  0  The  cutting  plane  intersects  the  THC  within  the  omni¬ 
directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between  a 
circular  arc  and  a  hyperbola. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  36.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  36  are  shown  in  Figure  37. 

First,  the  coordinates  of  the  ‘first  point’  on  the  THC,  p±,  are  determined  by 

_/  0  :7  +  M  =  f 

P1  \  tan  (7+|*|)  :  7  +  101  7^  |  ’  (44) 

Vpi  =  0. 
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(a)  |0|  +  7  <  tt/2 
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►  X 


(b)  |0|  +  7  >  tt/2 


Figure  35:  Case  A2A 


Similarly,  the  ‘first  point’  on  the  THC,  q\  is  determined  by 

:  7-  101  =  f 


i- <ji 


0 


v  tan (7-|0|)  :  "1  1^1  ^  2 

Vqi  =  0. 

It  can  be  shown  by  trigonometry  then  that  using  the  intermediate  value  g,  found  by  the  relation 

9=(b-v/^  +  <)^X, 

the  coordinates  of  the  ‘second  point,’  P2+-,  and  the  intermediate  value  q  can  be  determined  by 

£.r_r 

\h\ 


(45) 


(46) 


UP2+  —  XPl  T  9l 

q  = 


tan  |())| 


sin  |0|, 


(47) 


q2. 


yP2+  =  t/m2  -  < 

It  then  follows  that  the  hyperbolic  parameters  a,  6,  and  xs>  may  be  determined  by 

Xpi  Xfj 


vq  i 


b  =  a\/  e2  —  1, 


Xo7  —  U. 


(48) 


—  (49) 

In  the  implementation,  the  hyperbolic  span  is  split  into  two  regions,  such  that  the  initial  point  of  the 
outer  polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of 
visualization  functions.  The  upper  half  of  the  hyperbola  is  thus  spanned  by  the  angles 


0»  =  TT, 
ifjf  =  arctan 


Up2+ 


VP  2+ 


-  {x's  +  ae)  J  ’ 


(50) 
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relative  to  the  focus,  /.  The  radius  relative  to  /  is  determined  using  Equation  42.  The  circular  arc,  relative 
to  the  origin  is  spanned  by 


i pi  =  arctan 


(  VP2+ 


'P  2+ 


ipf  =  2tt  —  arctan 


f  VP2+ 


(51) 


\  XP2+  . 

with  a  radius  given  by  Equation  43.  Finally,  the  lower  half  of  the  hyperbola  is  spanned  by 

Up2+ 


if>i  =  —  arctan 
ij)f  =  —7 r. 


VP  2+ 


~  (x’s  +ae)J  : 


(52) 


Case  A2C:  0  >  h  and  h  >  — i?sin (7  +  |(^|)  The  cutting  plane  intersects  the  THC  within  the  omni¬ 
directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between  a 
circular  arc  and  a  hyperbola. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  38.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  38  are  shown  in  Figure  39. 

First,  the  coordinates  of  the  ‘first  point’  on  the  TF1C,  pi,  are  determined  by 

_/  0  : 7-  \<j>\  =  % 

Pl  \  tan  (7— 1</>|)  :  7  -  \4>\  7^  f  ’  (53) 

Upi  =  0. 

Similarly,  the  ‘first  point’  on  the  THC,  q\  is  determined  by 

I  0  : 7 +1^1  =  1 

qi  \  _tan(7+|0|)  :7+l<(,l^f  ’  (54) 

Vqi  =  0- 
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Figure  37:  Case  A2B 


It  can  be  shown  by  trigonometry  then  that  using  the  intermediate  values  q  and  g ,  found  by  the  relations 


9  = 

9  = 


(55) 


the  coordinates  of  the  ‘second  point,’  P2+5  can  be  determined  by 


XP2+  —  XPl  9 • 

yP2+  =  \]m2  -  q2. 

It  then  follows  that  the  hyperbolic  parameters  a,  6,  and  xsi  may  be  determined  by 

%pi  X q-i 

2  ’ 
b  =  a\/  e2  —  1, 


(56) 


(57) 


xs>  =  xPl  —  a.  (58) 

In  the  implementation,  the  hyperbolic  span  is  split  into  two  regions,  such  that  the  initial  point  of  the 
outer  polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of 
visualization  functions.  The  upper  half  of  the  hyperbola  is  thus  spanned  by  the  angles 


4>i  =  7T, 


ipf  =  arctan 


VP2+ _ 

-  K  +  ae) 


(59) 


relative  to  the  focus,  /.  The  radius  relative  to  /  is  determined  using  Equation  42.  The  circular  arc,  relative 
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to  the  origin  is  spanned  by 


(VP2+ 


'P2+ 


ij)i  =  arctan 
ipf  =  2n  —  arctan 

with  a  radius  given  by  Equation  43.  Finally,  the  lower  half  of  the  hyperbola  is  spanned  by 


f  VP2+ 

V  XP2+ 


ipi  =  —  arctan 
ifif  =  —7 r. 


VP2+ _ 

-  (x's  +  ae) 


(60) 


(61) 


z  z 


Case  A2D:  — i?sin(7+  \(f>\)  >  h  and  h  >  —R  The  cutting  plane  lies  below  the  lower  intersection  of 
the  THC  with  the  omni-directional  sensor  radius,  leading  to  a  disk-shaped  sensor  cross-section  with  radius 

^circular  arc- 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  40.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  40  are  shown  in  Figure  41. 

6.1.5  Case  Bl:  1  +  ne  >  e  >  1  —  ne  and  cf>  >  0 

The  THC  cross-section  is  strictly  parabolic,  thus,  sensor  region  cross-sections  may  be  disks,  non-existent, 
or  defined  by  a  compound  boundary  between  a  circular  arc  and  a  parabola.  Excluding  the  cases  where 
the  cross-section  is  determined  to  be  a  disk,  the  problem  is  determining  the  parabolic  parameter  P  in  the 
canonical  form  of  a  parabola,  and  the  ‘first  point,’  p\1  of  the  THC  conic  section.  ‘First  point’  refers  to  the 
point  on  the  parabola  most  negative  in  the  x  direction.  By  finding  the  location  of  a  ‘second  point,’  P2+, 
where  the  conic  section  intersects  the  omni-directional  sensor  radius  in  the  cutting  plane,  the  appropriate 
angular  intervals  in  ip  may  be  determined  to  define  the  hyperbolic  and  circular  curves. 
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(a)  |0|  +7  <  tt/2 


(b)  |0|  +  7  >  7r/2 


Figure  39:  Case  A2C 
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Figure  41:  Case  A2D 


For  the  parabolic  segments,  the  radius  of  the  parabola  with  respect  to  the  focus,  /,  located  on  the  x  axis 
at  xPl  +  P,  is  given  by 

2  p 

^parabola  =  Z  7  (62) 

1  —  COS  ip 

The  radius  of  the  outer  circular  arc,  bounding  the  omni-directional  sensor  radius  in  the  ATH  regions  is 
given  with  respect  to  the  origin  by 

X circular  arc  =  \J  R?  ~  h2 .  (63) 

Case  B1A:  R  >  h  and  h  >  0  The  cutting  plane  lies  above  the  upper  intersection  of  the  THC  with  the 
omni-directional  sensor  radius,  leading  to  a  disk-shaped  sensor  cross-section  with  radius  rc ircuiar  arc- 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  42.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  42  are  shown  in  Figure  43. 

Case  BIB:  0  >  h  and  h  >  —R sin27  The  cutting  plane  intersects  the  THC  within  the  omni-directional 
sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between  a  circular  arc  and 
a  parabola. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  44.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  44  are  shown  in  Figure  45. 

First,  the  coordinates  of  the  ‘first  point’  on  the  THC,  pi,  are  determined  by 

:27=f 

:  27  ^  |  ’  (64) 

Vpi  =  0- 


°Pl 


0 

\h\ 

tan  27 
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(a)  |0|  +  7  <  tt/2  (b)  |0|  +  7  >  tt/2 


Figure  43:  Case  B1A 
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It  can  be  shown  by  trigonometry  that  the  coordinates  of  the  ’second  point,’  P2+  can  be  determined  by 


VP2+ 


=  xPl  +  R  —  sqrthr  +  x: 


pi 


q  = 


-  (R  -  sqrth 2  +  x2Vl) 


sin  27 


cos  7 


(65) 


VP2+ 


=  \Jm2  —  q2. 


It  then  follows  that  the  parabolic  parameter  P  may  be  determined  by 

1,2 


P  = 


y\ 


p  2+ 


4  (xP2+  -  xpiy 


(66) 


In  the  implementation,  the  parabolic  span  is  split  into  two  regions,  such  that  the  initial  point  of  the  outer 
polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of  visual¬ 
ization  functions.  The  upper  half  of  the  parabola  is  thus  spanned  by  the  angles 


=  7T, 

ipf  =  arctan 


Vp2+ 


bp  2+ 


—  ( xPl  +  P)  J 


(67) 


relative  to  the  focus,  /.  The  radius  relative  to  /  is  determined  using  Equation  62.  The  circular  arc,  relative 
to  the  origin  is  spanned  by 


1 pi  =  arctan 


f  VP2+ 


'P  2+ 


ipf  =  2tt  —  arctan 


f  yP  2+ 
V  XP2+ 


(68) 


with  a  radius  given  by  Equation  63.  Finally,  the  lower  half  of  the  parabola  is  spanned  by 

VP2+ 


ipi  =  —  arctan 

tjjf  =  —  7T. 


up  2+ 


-  (Xpi  +P)J  ’ 


(69) 


Case  B1C:  — I?sin27  >  h  and  h  >  —R  The  cutting  plane  lies  below  the  lower  intersection  of  the  THC 
with  the  omni-directional  sensor  radius  leading  to  either  a  disk-shaped  sensor  cross-section  with  radius 
^ circular  arc,  or  no  region  at  all  if  101  +  7  >  7t/2. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  46.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  46  are  shown  in  Figure  47. 

6.1.6  Case  B2:  l+ne>e>  1  —  ne  and  cf>  <  0 

The  THC  cross-section  is  strictly  parabolic,  thus,  sensor  region  cross-sections  may  be  disks,  non-existent, 
or  defined  by  a  compound  boundary  between  a  circular  arc  and  a  parabola.  Excluding  the  cases  where 
the  cross-section  is  determined  to  be  a  disk,  the  problem  is  determining  the  parabolic  parameter  P  in  the 
canonical  form  of  a  parabola,  and  the  ‘first  point,’  p\ ,  of  the  THC  conic  section.  ‘First  point’  refers  to  the 
point  on  the  parabola  most  negative  in  the  x  direction.  By  finding  the  location  of  a  ‘second  point,’  P2+, 
where  the  conic  section  intersects  the  omni-directional  sensor  radius  in  the  cutting  plane,  the  appropriate 
angular  intervals  in  ip  may  be  determined  to  define  the  hyperbolic  and  circular  curves. 

For  the  parabolic  segments,  the  radius  of  the  parabola  with  respect  to  the  focus,  /,  located  on  the  x  axis 
at  xPl  +  P,  is  given  by 

2  p 

^parabola  —  Z  7  (70) 

1  —  COS  ip 

The  radius  of  the  outer  circular  arc,  bounding  the  omni-directional  sensor  radius  in  the  ATH  regions  is 
given  with  respect  to  the  origin  by 

^circular  arc  =  V  b" .  (71) 
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(a)  |0|  +7  <  V2 


(b)  |0|  +  7  >  7r/2 


Figure  44:  Case  BIB 
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Figure  45:  Case  BIB 
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(a)  |0|  +7  <  V2 


(b)  |0|  +  7  >  7r/2 


Figure  46:  Case  B1C 
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(b)  |0|  +  7  >  tt/2 


Figure  47:  Case  B1C 
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Case  B2A:  R  >  h  and  h  >=  R  sin  2y  The  cutting  plane  lies  above  the  upper  intersection  of  the  THC  with 
the  omni-directional  sensor  radius  leading  to  either  a  disk-shaped  sensor  cross-section  with  radius  rcircuial-  arc, 
or  no  region  at  all  if  \<j>\  +  7  >  7t/2. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  48.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  48  are  shown  in  Figure  49. 


(a)  |0|  +7  <  V2 


(b)  |0|  +  7  >  7I-/2 


Figure  48:  Case  B2A 


Case  B2B:  R  sin  27  >  h  and  h  >  0  The  cutting  plane  intersects  the  THC  within  the  omni-directional 
sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between  a  circular  arc  and 
a  parabola. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  50.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  50  are  shown  in  Figure  51. 

First,  the  coordinates  of  the  ‘first  point’  on  the  THC,  pi,  are  determined  by 


xv,  — 


0 

\h\ 


:27=f 

tan  27  :  27  7^  f 


(72) 


Vpi  =  0. 

It  can  be  shown  by  trigonometry  that  the  coordinates  of  the  ’second  point,’  P2+  can  be  determined  by 

xP2+  =  xPl  +  R  -  sqrth 2  +  xl  , 


<7  = 

yP2+  =  Vm2- 


,  „  ,0  9  ,  sin  2y 

m  —  [R  —  sqrth  +  xpi ) - 


cosy 


(73) 


q2. 


It  then  follows  that  the  parabolic  parameter  P  may  be  determined  by 

,,2 


P  = 


V; 


P  2+ 


4  (xP2+  -  xPl )' 


(74) 
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(a)  |0|  +  7  <  tt/2 


(b)  |0|  +  7  >  7r/2 


Figure  49:  Case  B2A 


In  the  implementation,  the  parabolic  span  is  split  into  two  regions,  such  that  the  initial  point  of  the  outer 
polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of  visual¬ 
ization  functions.  The  upper  half  of  the  parabola  is  thus  spanned  by  the  angles 


4>i  =  7T, 


ipf  =  arctan 


VP2+ _ 

—  ( xPl  +  P ) 


(75) 


relative  to  the  focus,  /.  The  radius  relative  to  /  is  determined  using  Equation  70.  The  circular  arc,  relative 
to  the  origin  is  spanned  by 


i pi  =  arctan 


VP2+ 

xP2+ 


ijif  =27 r  —  arctan 


with  a  radius  given  by  Equation  71.  Finally,  the  lower  half  of  the  parabola  is  spanned  by 


(76) 


ipi  =  —  arctan 
tf)f  =  —tt 


VP2+ _ 

—  ixpi  +  P) 


(77) 


Case  B2C:  0  >  h  and  h  >  —R  The  cutting  plane  lies  below  the  lower  intersection  of  the  THC  with  the 
omni-directional  sensor  radius,  leading  to  a  disk-shaped  sensor  cross-section  with  radius  rcjrcuiar  arc. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  52.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  52  are  shown  in  Figure  53. 
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(a)  |0|  +  7  <  V2 


(b)  |0|  +  7  >  7r/2 


Case  B2B 


(b)  |0|  +7  >  tt/2 


Figure  51:  Case  B2B 
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(a)  |0|  +  7  <  V2 


(b)  |0|  +  7  >  7r/2 


Figure  53:  Case  B2C 
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6.1.7  Case  Cl:  1  —  ne  >  e  and  <fi  >  0 


The  THC  cross-section  is  strictly  elliptical,  thus,  sensor  region  cross-sections  may  be  disks,  non-existent,  or 
defined  by  a  compound  boundary  between  a  circular  arc  and  an  ellipse  or  elliptical  arc.  Excluding  the  cases 
where  the  cross-section  is  determined  to  be  a  disk,  the  problem  is  determining  the  elliptical  parameters  a 
and  b  in  the  canonical  form  of  an  ellipse,  and  determining  the  center-point  of  the  conic-section,  s' .  s'  is 
centered  between  the  ‘first  point,’  pi,  of  the  THC  conic  section,  and  the  ‘third  point,’  P3,  at  the  opposite 
end  of  the  ellipse  major  axis  from  p\.  When  it  exists,  by  finding  the  location  of  a  ‘second  point,’  P2+,  where 
the  conic  section  intersects  the  omni-directional  sensor  radius  in  the  cutting  plane,  the  appropriate  angular 
intervals  in  if)  may  be  determined  to  define  the  hyperbolic  and  circular  curves.  When  p2+  does  not  exist,  the 
ellipse  forms  a  hole  in  a  disk  shaped  region. 

For  the  elliptical  segments,  the  radius  of  the  ellipse  with  respect  to  the  center-point  s' ,  is  given  by 

ab  . 

^ellipse  =  j  (78) 

J  ( b  cos  V’)2  +  (a  sin  i/j) 

The  radius  of  the  outer  circular  arc,  bounding  the  omni-directional  sensor  radius  in  the  ATH  regions  is 
given  with  respect  to  the  origin  by 

^circular  arc  =  R2  ~  h2 .  (79) 

Case  CIA:  R  >  h  and  h  >  0  The  cutting  plane  lies  above  the  upper  intersection  of  the  THC  with  the 
omni-directional  sensor  radius,  leading  to  a  disk-shaped  sensor  cross-section  with  radius  ?'circuiar  arc- 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  54.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  54  are  shown  in  Figure  55. 


Case  C1B:  0  >  h  and  h  >  —  i?sin(|</>|  —7)  The  cutting  plane  intersects  the  THC  within  the  omni¬ 
directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  circular  outer  boundary  with  an 
elliptical  interior  hole. 
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(a)  |0|  +  7  <  7r/2  (b)  |0|  +  7  >  7t/2 


Figure  55:  Case  CIA 


Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  56.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  56  are  shown  in  Figure  57. 

First,  the  coordinates  of  the  ‘first’  and  ‘third’  points,  p\  and  p$,  which  are  on  either  end  of  the  major 
axis  of  the  elliptical  THC  conic  section,  are  determined  by 

in 

te»(W+T>'  (80) 

\h\ 

tan  (|0|  +7)' 

It  then  follows  that  the  elliptical  parameters  a,  6,  and  x's  may  be  determined  by 


xPi  — 
Xp3  = 


Xp3  Xp^ 

2 

(81) 

a\J\  —  e2, 

Xpi  +  Xp3 

2 

(82) 

For  a  radius  given  by  equation  79,  an  outer  circle  spans  0  <  0  <  27 r,  while  an  internal  elliptical  hole  region 
is  defined  by  78  over  0  <i/j  <2tt  centered  on  s'. 


Case  C1C:  —R sin  (|0|  —  7)  >  h  and  h  >  — i?sin  (|0|  +  7)  The  cutting  plane  intersects  the  THC  within  the 
omni-directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between 
a  circular  arc  and  an  elliptical  arc. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  58.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  58  are  shown  in  Figure  59. 

First,  the  coordinates  of  the  ‘first’  and  ‘third’  points,  p\  and  p 3,  which  are  on  either  end  of  the  major 
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axis  of  the  elliptical  THC  conic  section,  are  determined  by 


\h\ 

tan  (|0|  +7)’ 

\h\ 

tan  (|0|  +  7) ' 


(83) 


It  can  be  shown  by  trigonometry  then  that  using  the  intermediate  values  q  and  g,  found  by  the  relations 

9  = 

<1  = 

cos  7 

the  coordinates  of  the  ‘second  point,’  P2+1  can  be  determined  by 


(R-y/ht  +  x^) 

sin  (7+  |0|) 


cos  7 
cos  \(j)\  ’ 


(84) 


Xp2_|_  —  %pi  H-  9 
VP2+  =  \/m2  -  q2 


(85) 


It  then  follows  that  the  elliptical  parameters  a,  b,  and  xsi  may  be  determined  by 


a  = 

*^Pl 

2 

(86) 

b  = 

a\J  1  —  e2, 

x's  = 

Xp1  ~b  Xp3 

(87) 

2 

In  the  implementation,  the  elliptical  arc  span  is  split  into  two  regions,  such  that  the  initial  point  of  the 
outer  polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of 
visualization  functions.  The  upper  half  of  the  elliptical  arc  is  thus  spanned  by  the  angles 


0i  =  7T, 


0/  =  arctan 


Up2+ 


Xp2+  X0 


(88) 


relative  to  s' .  The  radius  relative  to  s'  is  determined  using  Equation  78.  The  circular  arc,  relative  to  the 
origin  is  spanned  by 


{  VP2+ 

\a;P2+ 


0i  =  arctan 
0/  =  2n  —  arctan 


(  VP2+ 


T2+ 


with  a  radius  given  by  Equation  79.  Finally,  the  lower  half  of  the  hyperbola  is  spanned  by 


(89) 


0^  —  —  arctan 

0/  =  —7 r 


VP2+ 


(90) 


Case  C1D:  — i?sin(|0|  +  7)  >  h  and  h  >  —R  The  cutting  plane  lies  below  the  lower  intersection  of  the 
THC  with  the  omni-directional  sensor  radius  leading  to  either  a  disk-shaped  sensor  cross-section  with  radius 
‘r circular  arc,  or  no  region  at  all  if  101  +  7  >  7t/2. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  60.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  60  are  shown  in  Figure  61. 
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Figure  59:  Case  C1C 
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(b)  \c/>\  +7  >  ix / 2 


Figure  61:  Case  C1D 
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6.1.8  Case  C2:  1  —  ne  >  e  and  <p  <  0 


The  THC  cross-section  is  strictly  elliptical,  thus,  sensor  region  cross-sections  may  be  disks,  non-existent,  or 
defined  by  a  compound  boundary  between  a  circular  arc  and  an  ellipse  or  elliptical  arc.  Excluding  the  cases 
where  the  cross-section  is  determined  to  be  a  disk,  the  problem  is  determining  the  elliptical  parameters  a 
and  b  in  the  canonical  form  of  an  ellipse,  and  determining  the  center-point  of  the  conic-section,  s' .  s'  is 
centered  between  the  ‘first  point,’  pi,  of  the  THC  conic  section,  and  the  ‘third  point,’  P3,  at  the  opposite 
end  of  the  ellipse  major  axis  from  p\.  When  it  exists,  by  finding  the  location  of  a  ‘second  point,’  P2+,  where 
the  conic  section  intersects  the  omni-directional  sensor  radius  in  the  cutting  plane,  the  appropriate  angular 
intervals  in  ip  may  be  determined  to  define  the  hyperbolic  and  circular  curves.  When  p2+  does  not  exist,  the 
ellipse  forms  a  hole  in  a  disk  shaped  region. 

For  the  elliptical  segments,  the  radius  of  the  ellipse  with  respect  to  the  center-point  s',  is  given  by 

ab 

^ellipse  —  j  (91) 

J  ( b  cos  ip)  +  (a  sin  ip) 

The  radius  of  the  outer  circular  arc,  bounding  the  omni-directional  sensor  radius  in  the  ATH  regions  is 
given  with  respect  to  the  origin  by 

^circular  arc  =  \/  h  ~ .  (92) 

Case  C2A:  R  >  h  and  h  >  i?sin(|</>|  +7)  The  cutting  plane  lies  above  the  upper  intersection  of  the 
THC  with  the  omni-directional  sensor  radius  leading  to  either  a  disk-shaped  sensor  cross-section  with  radius 
t circular  arc,  or  no  region  at  all  if  |</>|  +  7  >  tt/2. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  62.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  62  are  shown  in  Figure  63. 


Case  C2B:  i? sin  ( |0|  +7 )  >  h  and  h  >  i?sin(|^|  —7)  The  cutting  plane  intersects  the  THC  within  the 
omni-directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  compound  boundary  between 
a  circular  arc  and  an  elliptical  arc. 
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(b)  |(/)|  +  7  >  tt/2 


Figure  63:  Case  C2A 


Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  64.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  64  are  shown  in  Figure  65. 

First,  the  coordinates  of  the  ‘first’  and  ‘third’  points,  p\  and  p$,  which  are  on  either  end  of  the  major 
axis  of  the  elliptical  THC  conic  section,  are  determined  by 


\h\ 


Xp 1  tan  (101 +7)’ 

a  -  H 


(93) 


P3  tan  (|0|  +  7)  ’ 

It  can  be  shown  by  trigonometry  then  that  using  the  intermediate  values  q  and  g7  found  by  the  relations 

9=(R~\/h2  +  x^)^W\’ 

sin  (7+  |0|) 


(94) 


q  =  m  —  g- 


COS7 


the  coordinates  of  the  ‘second  point,’  can  be  determined  by 


°p  2+ 


=  xPl  +  g. 


yP2+  =  V™2  ~  Q2- 

It  then  follows  that  the  elliptical  parameters  a,  b,  and  xs>  may  be  determined  by 

OCps 

O'  —  ~  1 


(95) 


b  =  a\J  1  — 


(96) 

(97) 
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In  the  implementation,  the  elliptical  arc  span  is  split  into  two  regions,  such  that  the  initial  point  of  the 
outer  polygon  boundary  is  at  the  satellite  location.  This  is  done  for  convenience  in  the  implementation  of 
visualization  functions.  The  upper  half  of  the  elliptical  arc  is  thus  spanned  by  the  angles 


Ipi  =  7T, 


4>f  =  arctan 


VP2+ 


Xp2+  Xa 


(98) 


relative  to  s' .  The  radius  relative  to  s'  is  determined  using  Equation  78.  The  circular  arc,  relative  to  the 
origin  is  spanned  by 


fVp2+ 

\XP2+ 


%l>i  =  arctan 
i/jf  =  2n  —  arctan 


(  VP2+ 


'P  2  + 


with  a  radius  given  by  Equation  79.  Finally,  the  lower  half  of  the  hyperbola  is  spanned  by 

,  |  VP2+ 

ipi  =  —  arctan  ' 
ipf  =  —7 r. 


(99) 


(100) 


z 
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Case  C2C:  iisin  (|0|  —  7)  >  h  and  h  0  The  cutting  plane  intersects  the  r  I "  1 1C'  wit  Inn  the  omm 
directional  sensor  radius,  leading  to  a  sensor  cross-section  defined  by  a  circular  outer  boundary  with  an 
elliptical  interior  hole. 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  66.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  66  are  shown  in  Figure  67. 
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(a)  |0|  +  7  <  tt/2 


y 


Figure  65:  Case  C2B 


First,  the  coordinates  of  the  ‘first’  and  ‘third’  points,  pi  and  p$,  which  are  on  either  end  of  the  major 
axis  of  the  elliptical  THC  conic  section,  are  determined  by 


=  \h\ 

tan  (|0|  +7)’ 

\h\ 

tan  (|0|  +7)' 


(101) 


It  then  follows  that  the  elliptical  parameters  a,  6,  and  x's  may  be  determined  by 


a  = 
b  = 


% P3  Xp  i 


(102) 


X 


/ 

s 


Xpi  H-  Xp3 
2 


(103) 


For  a  radius  given  by  equation  92,  an  outer  circle  spans  0  <  0  <  2-7t,  while  an  internal  elliptical  hole  region 
is  defined  by  91  over  0  <  ip  <2tt  centered  on  s'. 


Case  C2D:  0  >  h  and  h  >  — R  The  cutting  plane  lies  below  the  lower  intersection  of  the  THC  with  the 
omni-directional  sensor  radius,  leading  to  a  disk-shaped  sensor  cross-section  with  radius  rcircuiar  arc- 

Cutting  plane  positions  leading  to  this  case  are  shown  in  Figure  68.  The  associated  cutting-plane  cross- 
sections  of  the  ATH  coverage  volumes  in  Figure  68  are  shown  in  Figure  69. 


6.1.9  Transformation  from  xyz  to  xyz 

The  cross-sections  are  derived  in  a  Cartesian  coordinate  system  (xyz)  such  that  the  analysis  plane  is  parallel 
to  the  x  —  y  plane  in  the  general  problem,  the  z  axis  is  parallel  to  the  z  axis  in  the  general  problem,  and  the 
x  axis  points  toward  the  Earth  collincar  with  the  projection  of  the  satellite  position  vector  into  the  x  —  y 


68 


2 


2 


Figure  66:  Case  C2C 
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Figure  69:  Case  C2D 
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plane.  Transformation  from  the  xyz  frame  to  the  xyz  frame  is  accomplished  by  a  single  rotation  about  the  z 
axis  by  6  +  n,  where  9  is  the  longitude  of  the  satellite  location  in  the  xy  plane  (i.e.  angle  from  the  positive  x 
axis),  followed  by  a  simple  translation  of  coordinates  from  (0,  0, 0)  to  rsat,  where  rsat  is  the  satellite  position 
vector. 
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6.2  Target  Cross-Section  Generation 

Target  regions  are  typically  geometrically  less  complex  than  an  ATH  sensor  region  which  consists  of  a 
sphere  and  a  missing  cone  region.  Thus  target  region  cross-sections  are  significantly  more  straightforward 
to  compute  in  most  cases. 

6.2.1  Dual- Altitude  Band  Target  Volume 

The  target  region  considered  in  the  current  study  is  a  simple  dual-altitude  band  target  volume  -  a  hollow 
sphere  with  an  upper  target  radius  of  ru,  and  a  lower  target  radius  of  n,  such  that  r*  <  77  <  ru.  Figure  70 
shows  several  cross-sections  of  a  dual-altitude  band  target  volume. 


Figure  70:  Offset  cross-sections  of  a  dual-altitude  band  target  volume 


Case  A:  \h\  <  ri  The  cutting  plane  intersects  shells  formed  by  both  the  outer  and  inner  target  radii,  thus 
the  cross-section  is  an  annulus  with  an  outer  radius 

ru>  =  Vrl~  h2,  (104) 

and  an  inner  radius 

7 v  =  \J  rf  —  h2 .  (105) 

Case  B:  n  <  \h\  <  ru  The  cutting  plane  intersects  only  the  shell  formed  by  the  outer  target  radius,  thus 
the  cross-section  is  a  disk  with  radius 

rv  =  Vrl  -  h2-  (106) 

Case  C:  ru  <  \h\  The  cutting  plane  does  not  intersect  the  target  region  at  all,  and  thus  no  cross-section 
exists  at  this  value  of  h. 

6.2.2  Other  Target  Volume  Types 

Any  target  volume  for  which  cross-section  generation  is  straightforward  may  be  implemented  in  an  analysis 
with  minimal  effort.  For  instance,  the  dual- altitude  band  target  volume  case  may  be  modified  to  consider 
lower  and  upper  latitude  limits.  The  resulting  target  volume  would  be  toroidal  in  shape.  Similarly,  an 
Earth-fixed  target  volume  may  be  defined  incorporating  latitude,  longitude,  and  altitude  bounds  to  describe 
a  volume  of  space  fixed  above  some  geographic  region. 
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7  Non-Planar  Analysis  Example  Problems 

The  Examples  presented  in  this  section  focus  on  analysis  only  at  this  time.  The  methods  developed  in  the 
previous  section  are  thus  far  only  implemented  in  MATLAB,  using  the  MATLAB/C  version  of  the  planar 
analysis  algorithm.  Thus,  performance  is  insufficient  for  the  optimization  of  any  practical  constellation  in  a 
reasonable  amount  of  computer  time  executing  on  a  single  desktop  computer.  This  issue  will  be  alleviated 
in  future  work  when  the  C/C++  implementation  is  parallelized  and  updated  to  analyze  the  volumetric  case. 

7.1  Example  7  —  ATH  Coverage  by  a  Small  Non-Coplanar  Constellation 

A  constellation  is  considered  with  the  properties  shown  in  Table  13.  The  parameter  hles  refers  to  the 
number  of  equally  spaced  offset  planes  in  which  cross-sectional  coverage  area  is  computed.  These  planes  are 
distributed  between  on  the  interval  h  €  [— (re  +  hu ),  (re  +  hu)\. 

Table  13:  Example  7  -  Parameters 


Parameter 

Description 

Value 

ht 

tangent  height 

100  km 

hi 

lower  target  altitude 

1000  km 

hu 

upper  target  altitude 

5000  km 

k 

number  of  orbital  planes 

3 

rii 

satellites  per  plane 

2 

cii 

semi- major  axis 

10000  km 

orbit  eccentricity 

0.25° 

ii 

inclination  of  plane  1 

0° 

*2 

inclination  of  plane  2 

45° 

h 

inclination  of  plane  3 

90° 

RAAN  of  plane  1 

0° 

02 

RAAN  of  plane  2 

120° 

fb 

RAAN  of  plane  3 

240° 

LOi 

argument  of  periapsis  of  each  orbit 

0° 

A  Mi 

spacing  in  mean  anomaly 

180° 

Mq  i 

mean  anomaly  of  lead  satellite  at  epoch 

0° 

R 

omni-directional  sensor  range 

5000  km 

m 

initial  polygon  resolution 

50  PPC 

^•res 

offset  plane  resolution 

100 

The  single  and  double  ATH  coverage  provided  by  this  constellation  over  3  periods  is  shown  in  Figure  71. 
See  accompanying  video  #4  for  an  animation  of  the  constellation  and  coverage  evolution. 

The  execution  time  of  this  analysis  (at  500  time  steps)  in  MATLAB  is  approximately  20  minutes  in  64-bit 
Windows  7  on  a  3.6GHz  Intel  Core  i7  platform. 

7.2  Example  8  —  ATH  Coverage  by  a  Walker  Star  Constellation 

The  Walker  Star  constellation  is  a  traditional  satellite  constellation  for  delivering  BTH  coverage.  The  Iridium 
satellite  phone  constellation  is  one  operational  example  of  a  Walker  Star.  Idealized  parameters  for  the  Iridium 
constellation  are  shown  in  Table  14.  The  parameter  hles  refers  to  the  number  of  equally  spaced  offset  planes 
in  which  cross-sectional  coverage  area  is  computed.  These  planes  are  distributed  between  on  the  interval 
h  (E  [  (re  +  hu),  (re  + 

The  single  and  double  ATH  coverage  provided  by  this  constellation  over  20  minutes  is  shown  in  Figure 
72.  See  accompanying  video  ^5  for  an  animation  of  the  constellation  and  coverage  evolution. 

The  origin  of  the  Walker  Star  constellation  is  in  the  ‘streets  of  coverage’  approach,  where  satellite  planes 
are  spaced  and  phased  such  that  continuous  single  BTH  coverage  is  available  at  the  equator  where  longitudinal 
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(b)  2x  Coverage 


Figure  71:  Example  7  Single  and  Double  ATH  Coverage  (see  accompanying  video  #4) 
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(a)  lx  Coverage 


Time  (sidereal  days) 


(b)  2x  Coverage 


Figure  72:  Example  8  -  Single  and  Double  ATH  Coverage  (see  accompanying  video  #5) 
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Table  14:  Example  8  -  Parameters 


Parameter 

Description 

Value 

ht 

tangent  height 

100  km 

hi 

lower  target  altitude 

200  km 

hu 

upper  target  altitude 

2500  km 

k 

number  of  orbital  planes 

6 

rii 

satellites  per  plane 

11 

hi 

satellite  altitude 

781  km 

e* 

orbit  eccentricity 

0° 

k 

plane  inclination 

86.4° 

fil 

RAAN  of  plane  1 

0° 

D‘2 

RAAN  of  plane  2 

31.6° 

C,3 

RAAN  of  plane  3 

63.2° 

^4 

RAAN  of  plane  4 

94.8° 

C.5 

RAAN  of  plane  5 

126.4° 

^6 

RAAN  of  plane  6 

158° 

A  Mi 

spacing  in  mean  anomaly 

32.73° 

PO 

phase  offset  b/w  adjacent  planes 

18° 

R 

omni-directional  sensor  range 

2500  km 

m 

initial  polygon  resolution 

50  PPC 

hres 

offset  plane  resolution 

100 

spacing  of  satellites  is  at  a  maximum.  From  this  example  it  is  clear  to  see  that  a  volumetric  analog  of  the 
‘streets  of  coverage’  approach  may  be  used  to  design  constellations  for  ATH  coverage  under  certain  conditions. 
The  definition  of  a  ‘street’  would  necessarily  be  amended  to  a  three-dimensional  analog  of  a  ‘street,’  such 
that  coverage  throughout  the  altitude  range  is  ensured  at  the  equator  in  the  final  solution. 

The  execution  time  of  this  analysis  (at  500  time  steps)  in  MATLAB  is  approximately  60  minutes  in  64-bit 
Windows  7  on  a  3.6GHz  Intel  Core  i7  platform. 

7.3  Example  9  —  ATH  Coverage  by  a  Walker  Delta  Constellation 

The  Walker  Delta  constellation  is  a  traditional  satellite  constellation  for  delivering  BTH  coverage.  The 
GlobalStar  satellite  phone  constellation  is  one  operational  example  of  a  Walker  Delta.  Idealized  parameters 
for  the  GlobalStar  constellation  are  shown  in  Table  15.  The  parameter  hles  refers  to  the  number  of  equally 
spaced  offset  planes  in  which  cross-sectional  coverage  area  is  computed.  These  planes  are  distributed  between 
on  the  interval  h  £  [— (re  +  hu),  (re  +  hu)\. 

The  single  and  double  ATH  coverage  provided  by  this  constellation  over  20  minutes  is  shown  in  Figure 
73.  See  accompanying  video  ^6  for  an  animation  of  the  constellation  and  coverage  evolution. 

Just  as  in  the  BTH  case,  this  low  altitude  Walker  Delta  constellation  provides  no  ATH  coverage  at  the 
poles.  Thus,  a  Walker  Delta  may  only  be  suitable  for  ATH  coverage  of  a  target  region  limited  to  some 
latitude  range  centered  at  the  equator. 

The  execution  time  of  this  analysis  (at  500  time  steps)  in  MATLAB  is  approximately  45  minutes  in  64-bit 
Windows  7  on  a  3.6GHz  Intel  Core  i7  platform. 
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(a)  lx  Coverage 


(b)  2x  Coverage 


Figure  73:  Example  9  Single  and  Double  ATH  Coverage  (see  accompanying  video  #6) 
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Table  15:  Example  9  Parameters 


Parameter 

Description 

Value 

ht 

tangent  height 

100  km 

hi 

lower  target  altitude 

200  km 

h"u 

upper  target  altitude 

2500  km 

k 

number  of  orbital  planes 

8 

rii 

satellites  per  plane 

6 

hi 

satellite  altitude 

1414  km 

e-i 

orbit  eccentricity 

0° 

k 

plane  inclination 

52.0° 

fil 

RAAN  of  plane  1 

0° 

RAAN  of  plane  2 

45° 

O3 

RAAN  of  plane  3 

90° 

RAAN  of  plane  4 

135° 

O5 

RAAN  of  plane  5 

180° 

06 

RAAN  of  plane  6 

225° 

RAAN  of  plane  7 

270° 

o8 

RAAN  of  plane  8 

315° 

AM, 

spacing  in  mean  anomaly 

60° 

PO 

phase  offset  b/w  adjacent  planes 

7.5° 

R 

omni-directional  sensor  range 

2500  km 

m 

initial  polygon  resolution 

50  PPC 

^res 

offset  plane  resolution 

100 
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8  Conclusions 


During  the  course  of  this  effort,  optimal  constellation  design  for  above-the- horizon  coverage  was  explored. 
Initially,  the  analysis  focused  only  on  planar  constellations.  This  facilitated  some  analytical  studies  that 
served  as  a  basis  for  validation  of  a  greater  numerical  process,  later  extended  to  address  a  simplified  version 
of  the  three-dimensional  problem.  The  work  yielded  two  significant  developments.  First,  a  generalized 
numerical  algorithm  was  designed,  tested,  and  validated  that  can  be  used  to  address  the  design  of  planar 
constellations  for  optimal  above-the-liorizon  coverage.  The  algorithm  allows  for  the  consideration  of  non- 
omnidirectional  sensors  and  arbitrary  sensor  pointing.  The  product  of  this  work  served  as  a  stepping  stone  to 
the  second  development,  the  extension  to  the  non-coplanar  (i.e.  three-dimensional)  case,  where  the  satellites 
in  the  constellation  exist  in  differing  elliptical  three-dimensional  orbits.  This  latest  aspect  of  the  development 
was  restricted  to  omni-directional  sensors,  for  simplicity,  but  the  overall  process  of  computing  the  cost  metric 
in  the  three-dimensional  case  was  otherwise  successfully  demonstrated. 

From  a  software  development  perspective,  both  the  two-  and  three-dimensional  cases  were  developed 
in  Matlab,  though  the  two-dimensional  case  also  explored  C/C++  as  a  means  of  increasing  computational 
efficiency.  Future  studies  may  consider  parallelization,  via  GPU’s,  of  the  overall  cost  metric  calculation, 
which  should  offer  significant  performance  improvements  in  any  optimization  process  that  depends  on  this 
cost  metric.  The  nature  of  the  problem  lends  itself  well  to  parallel  computing.  With  a  high-performance 
computing  implementation,  various  NLP  and  MINLP  solvers  may  be  wrapped  around  the  new  numerical 
volumetric  ATH  coverage  model.  Optimal  constellation  design  may  then  be  approached  in  a  volumetric 
sense  as  readily  as  it  was  addressed  in  the  planar  case  during  the  initial  years  of  this  award. 
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